WIP: assignment 2

This commit is contained in:
2024-11-25 01:21:43 +02:00
parent 91605e6f3a
commit 9f8702729d
15 changed files with 649 additions and 0 deletions
+22
View File
@@ -0,0 +1,22 @@
% Given environment
clear;
% Setup the function under test
syms x y;
fexpr = x^5 * exp(-x^2 - y^2);
title_fun = "$f(x,y) = x^5 \cdot e^{-x^2 - y^2}$";
% Calculate the gradient and Hessian
grad_fexpr = gradient(fexpr, [x, y]); % Gradient of f
hessian_fexpr = hessian(fexpr, [x, y]); % Hessian of f
% Convert symbolic expressions to MATLAB functions
fun = matlabFunction(fexpr, 'Vars', [x, y]); % Function
grad_fun = matlabFunction(grad_fexpr, 'Vars', [x, y]); % Gradient
hessian_fun = matlabFunction(hessian_fexpr, 'Vars', [x, y]); % Hessian
% Amijo globals
global amijo_beta amijo_sigma
%fixed step size globals
global gamma_fixed_step
+54
View File
@@ -0,0 +1,54 @@
% Define environment (functions, gradients etc...)
GivenEnv
% Define parameters
max_iter = 300; % Maximum iterations
tol = 1e-4; % Tolerance
% Methods tuning
amijo_beta = 0.5; % Armijo reduction factor
amijo_sigma = 0.1; % Armijo condition constant
m = 0.01;
% Point x0 = (0, 0)
% x0 = [0, 0];
% Point x0 = (0, 0)
x0s = [0, 0; -1, 1 ; 1, -1]; % Initial points
for i = 1:size(x0s, 1)
x0 = x0s(i, :);
point_str = "[" + x0(1) + ", " + x0(2) + "]";
% Find the best fixed gamma
k = zeros(100, 1);
j = 1;
n = linspace(0.1, 1.5, 100);
for g = n
gamma_fixed_step = g;
[~, ~, k(j)] = lev_mar(fun, grad_fun, hessian_fun, m, x0, tol, max_iter, 'fixed');
j = j + 1;
end
plotIterationsOverGamma(n, k, "Iteration for different $\gamma$ values", "LevMar_Iter_o_gamma_" + i + ".png");
[~, j] = min(k);
gamma_fixed_step = n(j);
[x_fixed, f_fixed, kk] = lev_mar(fun, grad_fun, hessian_fun, m, x0, tol, max_iter, 'fixed');
fprintf('Fixed step: Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_fixed(end, :), f_fixed(end));
plotPointsOverContour(x_fixed, fun, [-2, 2], [-3, 3], 100, point_str + ": LevMar $\gamma$ = " + gamma_fixed_step, "LevMar_fixed_" + i + ".png");
% Minimized f
[x_minimized, f_minimized, kk] = lev_mar(fun, grad_fun, hessian_fun, m, x0, tol, max_iter, 'minimized');
fprintf('Minimized f(g): Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_minimized(end, :), f_minimized(end));
plotPointsOverContour(x_minimized, fun, [-2, 2], [-3, 3], 100, point_str + ": LevMar minimized $f(x_k + \gamma_kd_k)$", "LevMar_minimized_" + i + ".png");
% Armijo Rule
[x_armijo, f_armijo, kk] = lev_mar(fun, grad_fun, hessian_fun, m, x0, tol, max_iter, 'armijo');
fprintf('Armijo step: Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_armijo(end, :), f_armijo(end));
plotPointsOverContour(x_armijo, fun, [-2, 2], [-3, 3], 100, point_str + ": LevMar Armijo method", "LevMar_armijo_" + i + ".png");
end
+52
View File
@@ -0,0 +1,52 @@
% Define environment (functions, gradients etc...)
GivenEnv
% Define parameters
max_iter = 300; % Maximum iterations
tol = 1e-4; % Tolerance
% Methods tuning
amijo_beta = 0.5; % Armijo reduction factor
amijo_sigma = 0.1; % Armijo condition constant
% Point x0 = (0, 0)
% x0 = [0, 0];
% Point x0 = (0, 0)
x0s = [0, 0; -1, 1 ; 1, -1]; % Initial points
for i = 1:size(x0s, 1)
x0 = x0s(i, :);
point_str = "[" + x0(1) + ", " + x0(2) + "]";
% Find the best fixed gamma
k = zeros(100, 1);
j = 1;
n = linspace(0.1, 1.5, 100);
for g = n
gamma_fixed_step = g;
[~, ~, k(j)] = newton(fun, grad_fun, hessian_fun, x0, tol, max_iter, 'fixed');
j = j + 1;
end
plotIterationsOverGamma(n, k, "Iteration for different $\gamma$ values", "Newton_Iter_o_gamma_" + i + ".png");
[~, j] = min(k);
gamma_fixed_step = n(j);
[x_fixed, f_fixed, kk] = newton(fun, grad_fun, hessian_fun, x0, tol, max_iter, 'fixed');
fprintf('Fixed step: Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_fixed(end, :), f_fixed(end));
plotPointsOverContour(x_fixed, fun, [-2, 2], [-3, 3], 100, point_str + ": Newton $\gamma$ = " + gamma_fixed_step, "Newton_fixed_" + i + ".png");
% Minimized f
[x_minimized, f_minimized, kk] = newton(fun, grad_fun, hessian_fun, x0, tol, max_iter, 'minimized');
fprintf('Minimized f(g): Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_minimized(end, :), f_minimized(end));
plotPointsOverContour(x_minimized, fun, [-2, 2], [-3, 3], 100, point_str + ": Newton minimized $f(x_k + \gamma_kd_k)$", "Newton_minimized_" + i + ".png");
% Armijo Rule
[x_armijo, f_armijo, kk] = newton(fun, grad_fun, hessian_fun, x0, tol, max_iter, 'armijo');
fprintf('Armijo step: Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_armijo(end, :), f_armijo(end));
plotPointsOverContour(x_armijo, fun, [-2, 2], [-3, 3], 100, point_str + ": Newton Armijo method", "Newton_armijo_" + i + ".png");
end
+52
View File
@@ -0,0 +1,52 @@
% Define environment (functions, gradients etc...)
GivenEnv
% Define parameters
max_iter = 300; % Maximum iterations
tol = 1e-4; % Tolerance
% Methods tuning
amijo_beta = 0.5; % Armijo reduction factor
amijo_sigma = 0.1; % Armijo condition constant
% Point x0 = (0, 0)
% x0 = [0, 0];
% Point x0 = (0, 0)
x0s = [-1, 1 ; 1, -1]; % Initial points
for i = 1:size(x0s, 1)
x0 = x0s(i, :);
point_str = "[" + x0(1) + ", " + x0(2) + "]";
% Find the best fixed gamma
k = zeros(100, 1);
j = 1;
n = linspace(0.1, 1.5, 100);
for g = n
gamma_fixed_step = g;
[~, ~, k(j)] = steepest_descent(fun, grad_fun, x0, tol, max_iter, 'fixed');
j = j + 1;
end
plotIterationsOverGamma(n, k, "Iteration for different $\gamma$ values", "StDes_Iter_o_gamma_" + i + ".png");
[~, j] = min(k);
gamma_fixed_step = n(j);
[x_fixed, f_fixed, kk] = steepest_descent(fun, grad_fun, x0, tol, max_iter, 'fixed');
fprintf('Fixed step: Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_fixed(end, :), f_fixed(end));
plotPointsOverContour(x_fixed, fun, [-2, 2], [-2, 2], 100, point_str + ": Steepest descent $\gamma$ = " + gamma_fixed_step, "StDes_fixed_" + i + ".png");
% Minimized f
[x_minimized, f_minimized, kk] = steepest_descent(fun, grad_fun, x0, tol, max_iter, 'minimized');
fprintf('Minimized f(g): Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_minimized(end, :), f_minimized(end));
plotPointsOverContour(x_minimized, fun, [-2, 2], [-2, 2], 100, point_str + ": Steepest descent minimized $f(x_k + \gamma_kd_k)$", "StDes_minimized_" + i + ".png");
% Armijo Rule
[x_armijo, f_armijo, kk] = steepest_descent(fun, grad_fun, x0, tol, max_iter, 'armijo');
fprintf('Armijo step: Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_armijo(end, :), f_armijo(end));
plotPointsOverContour(x_armijo, fun, [-2, 2], [-2, 2], 100, point_str + ": Steepest descent Armijo method", "StDes_armijo_" + i + ".png");
end
+26
View File
@@ -0,0 +1,26 @@
function [gamma] = gamma_armijo(f, grad_f, x0)
% Calculates the best step based on amijo method
%
% f(xk γk*f(xk)) f(xk) σ*γk*f(xk)^2
%
% f: Objective function
% x0: Initial (x,y) point
% beta: beta factor in (0, 1)
% signam: sigma factor in (0,1)
global amijo_beta
global amijo_sigma
gamma = 1; % Start with a step size of 1
grad = grad_f(x0(1), x0(2));
% Perform Armijo line search
while f(x0(1) - gamma * grad(1), x0(2) - gamma * grad(2)) > ...
f(x0(1), x0(2)) - amijo_sigma * gamma * norm(grad)^2
gamma = amijo_beta * gamma; % Reduce step size
end
end
+12
View File
@@ -0,0 +1,12 @@
function [gamma] = gamma_fixed(~, ~, ~)
% Return a fixed step
%
% This is for completion and code symmetry.
%
global gamma_fixed_step
% Perform line search
gamma = gamma_fixed_step;
end
+19
View File
@@ -0,0 +1,19 @@
function [gamma] = gamma_minimized(f, grad_f, x0)
% Calculates the step based on minimizing f(xk γ*f(xk))
%
%
% f: Objective function
% grad_f: Gradient of objective function
% x0: Initial (x,y) point
% Define the line search function g(gamma) = f(x0 - gamma * grad)
grad = grad_f(x0(1), x0(2));
g = @(gamma) f(x0(1) - gamma * grad(1), x0(2) - gamma * grad(2));
% Perform line search
gamma = fminbnd(g, 0, 1);
% ToDo: Check if we can use fmin_bisection_der
% from the previous assigment here!
end
+49
View File
@@ -0,0 +1,49 @@
function [x_vals, f_vals, k] = lev_mar(f, grad_f, hessian_f, m, x0, tol, max_iter, mode)
% f: Objective function
% grad_f: Gradient of the function
% hessian_f: Hessian of the function
% x0: Initial point [x0, y0]
% tol: Tolerance for stopping criterion
% max_iter: Maximum number of iterations
% x_vals: Vector with the (x,y) values until minimum
% f_vals: Vector with f(x,y) values until minimum
% k: Number of iterations
if strcmp(mode, 'armijo') == 1
gamma_f = @(f, grad_f, x0) gamma_armijo(f, grad_f, x0);
elseif strcmp(mode, 'minimized') == 1
gamma_f = @(f, grad_f, x0) gamma_minimized(f, grad_f, x0);
else % mode == 'fixed'
gamma_f = @(f, grad_f, x0) gamma_fixed(f, grad_f, x0);
end
x_vals = x0; % Store iterations
f_vals = f(x0(1), x0(2));
for k = 1:max_iter
grad = grad_f(x0(1), x0(2));
% Check for convergence
if norm(grad) < tol
break;
end
hess = hessian_f(x0(1), x0(2));
mI = m * eye(size(hess));
% Solve for search direction using Newton's step
dk = - inv(hess + mI) * grad;
% Calculate gamma
gamma = gamma_f(f, grad_f, x0);
x_next = x0 + gamma * dk'; % Update step
f_next = f(x_next(1), x_next(2));
x0 = x_next; % Update point
x_vals = [x_vals; x_next]; % Store values
f_vals = [f_vals; f_next]; % Store function values
end
end
+48
View File
@@ -0,0 +1,48 @@
function [x_vals, f_vals, k] = newton(f, grad_f, hessian_f, x0, tol, max_iter, mode)
% f: Objective function
% grad_f: Gradient of the function
% hessian_f: Hessian of the function
% x0: Initial point [x0, y0]
% tol: Tolerance for stopping criterion
% max_iter: Maximum number of iterations
% x_vals: Vector with the (x,y) values until minimum
% f_vals: Vector with f(x,y) values until minimum
% k: Number of iterations
if strcmp(mode, 'armijo') == 1
gamma_f = @(f, grad_f, x0) gamma_armijo(f, grad_f, x0);
elseif strcmp(mode, 'minimized') == 1
gamma_f = @(f, grad_f, x0) gamma_minimized(f, grad_f, x0);
else % mode == 'fixed'
gamma_f = @(f, grad_f, x0) gamma_fixed(f, grad_f, x0);
end
x_vals = x0; % Store iterations
f_vals = f(x0(1), x0(2));
for k = 1:max_iter
grad = grad_f(x0(1), x0(2));
% Check for convergence
if norm(grad) < tol
break;
end
hess = hessian_f(x0(1), x0(2));
% Solve for search direction using Newton's step
dk = - inv(hess) * grad;
% Calculate gamma
gamma = gamma_f(f, grad_f, x0);
x_next = x0 + gamma * dk'; % Update step
f_next = f(x_next(1), x_next(2));
x0 = x_next; % Update point
x_vals = [x_vals; x_next]; % Store values
f_vals = [f_vals; f_next]; % Store function values
end
end
+42
View File
@@ -0,0 +1,42 @@
function plot3Dfun(fun, x_lim, y_lim, size, plot_title)
% 3D plots a function
% fun: The function to plot
% x_lim: The range for x axis. ex: [-2, 2]
% y_lim: The range for y axis. ex: [0, 2]
% size: The number of points for each axis
% plot_title: The latex title for the plot
%
% Generate a grid for x and y
x_space = linspace(x_lim(1), x_lim(2), size);
y_space = linspace(y_lim(1), y_lim(2), size);
[X, Y] = meshgrid(x_space, y_space);
% Evaluate the function on the grid
Z = fun(X, Y);
% 3D plot
figure('Name', 'f(x,y)', 'NumberTitle', 'off');
set(gcf, 'Position', [100, 100, 960, 960]); % Set the figure size
surf(X, Y, Z);
% Customize the plot
xlabel('x'); % Label for x-axis
ylabel('y'); % Label for y-axis
zlabel('f(x, y)'); % Label for z-axis
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
colorbar;
% save the figure
print(gcf, 'FunctionPlot.png', '-dpng', '-r300');
% Contours
figure('Name', 'Contours of f(x,y)', 'NumberTitle', 'off');
set(gcf, 'Position', [100, 100, 1280, 1280]); % Set the figure size
contour(X, Y, Z);
xlabel('x'); % Label for x-axis
ylabel('y'); % Label for y-axis
title(plot_title, 'Interpreter', 'latex', 'FontSize', 20); % Title of the plot
colorbar;
print(gcf, 'ContoursPlot.png', '-dpng', '-r300');
end
+19
View File
@@ -0,0 +1,19 @@
function plotIterationsOverGamma(gamma, iterations, plot_title, filename)
% 3D plots a function
% fun: The points to plot
% contur_fun: The function for contour plot
% x_lim: The range for x axis. ex: [-2, 2]
% y_lim: The range for y axis. ex: [0, 2]
% size: The number of points for each axis
% plot_title: The latex title for the plot
% filename: The filename to save the plot (if exists)
%
figure('Name', 'Iterations_over_gamma', 'NumberTitle', 'off');
set(gcf, 'Position', [100, 100, 960, 960]); % Set the figure size
plot(gamma, iterations, '*r', 'LineWidth', 2);
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
xlabel('\gamma') ;
ylabel('Iterations');
print(gcf, filename, '-dpng', '-r300');
end
+40
View File
@@ -0,0 +1,40 @@
function plotPointsOverContour(points, contour_fun, x_lim, y_lim, size, plot_title, filename)
% 3D plots a function
% points: The points to plot
% contur_fun: The function for contour plot
% x_lim: The range for x axis. ex: [-2, 2]
% y_lim: The range for y axis. ex: [0, 2]
% size: The number of points for each axis
% plot_title: The latex title for the plot
% filename: The filename to save the plot (if exists)
%
% Generate a grid for x and y
x_space = linspace(x_lim(1), x_lim(2), size);
y_space = linspace(y_lim(1), y_lim(2), size);
[X, Y] = meshgrid(x_space, y_space);
% Evaluate the function on the grid
Z = contour_fun(X, Y);
% 2D plot
figure('Name', '(x,y)', 'NumberTitle', 'off');
set(gcf, 'Position', [100, 100, 960, 960]); % Set the figure size
plot(points(:, 1), points(:, 2), '-or');
hold on
contour(X, Y, Z);
% Customize the plot
xlim(x_lim);
ylim(y_lim);
xlabel('x'); % Label for x-axis
ylabel('y'); % Label for y-axis
grid on
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
colorbar;
% save the figure
if strcmp(filename, '') == 0
print(gcf, filename, '-dpng', '-r300');
end
end
+44
View File
@@ -0,0 +1,44 @@
function [x_vals, f_vals, k] = steepest_descent(f, grad_f, x0, tol, max_iter, mode)
% f: Objective function
% grad_f: Gradient of the function
% x0: Initial point [x0, y0]
% tol: Tolerance for stopping criterion
% max_iter: Maximum number of iterations
% x_vals: Vector with the (x,y) values until minimum
% f_vals: Vector with f(x,y) values until minimum
% k: Number of iterations
if strcmp(mode, 'armijo') == 1
gamma_f = @(f, grad_f, x0) gamma_armijo(f, grad_f, x0);
elseif strcmp(mode, 'minimized') == 1
gamma_f = @(f, grad_f, x0) gamma_minimized(f, grad_f, x0);
else % mode == 'fixed'
gamma_f = @(f, grad_f, x0) gamma_fixed(f, grad_f, x0);
end
% Storage for iterations, begin with the first point
x_vals = x0;
f_vals = f(x0(1), x0(2));
for k = 1:max_iter
grad = grad_f(x0(1), x0(2));
% Check for convergence
if norm(grad) < tol
break;
end
dk = - grad;
% Calculate gamma
gk = gamma_f(f, grad_f, x0);
x_next = x0 + gk * dk'; % Update step
f_next = f(x_next(1), x_next(2));
x0 = x_next; % Update point
x_vals = [x_vals; x_next]; % Store values
f_vals = [f_vals; f_next]; % Store function values
end
end