Question

To translate the steepest ascent into matlab code:

f(x,y)=-8*x + 12*y + x^2 - 2*x^4 - 2*x*y + 4*y^2

dx=-8+2x-8x^3-2y

dy=12-2x+8y

subs in x0=0, y0=0 into dx and dy

dx=-8

dy=12

subs dx and dy into x1=x0 + dx*0.001 and y1=y0 + dy*0.001

then subs x1 and y1 back to x0 and x0; repeat the step until the answer has reach optimisation.

Answer 1

clc
clear
close all

f= @(x,y) -8.*x + 12.*y + x.^2 - 2.*x.^4 - 2.*x.*y + 4.*y.^2;
dx = @(x,y) -8+2*x-8*x^3-2*y;
dy = @(x,y) 12-2*x+8*y;

x(1) = -4; y(1) = -4; %start at [-4,-4]
z(1) = f(x(1),y(1));

max_iter = 300;
for i=1:max_iter

x(i+1) = x(i) + 0.001*dx(x(i),y(i));
y(i+1) = y(i) + 0.001*dy(x(i),y(i));

z(i+1) = f(x(i+1),y(i+1));
diff(i) = z(i+1) - z(i);

end

optimum_iteration = find(diff==min(diff))
optimum_x = x(optimum_iteration)
optimum_y = y(optimum_iteration)
optimum_f = f(optimum_x,optimum_y)

figure
plot3(x,y,z)
hold on
plot3(optimum_x,optimum_y,optimum_f,'x r')
title('f(x,y) = -8*x + 12*y + x^2 - 2*x^4 - 2*x*y + 4*y^2')
xlabel('x');
ylabel('y')
zlabel('f(x,y)')
grid on
legend('f(x,y)','optimum f(x,y)')

figure
plot(diff,'linewidth',2)
hold on
plot(optimum_iteration,diff(optimum_iteration),'x r')
title('Change in f(x,y) with iteration')
ylabel('f(iter+1) - f(iter)')
xlabel('iteration (iter)')
legend('change in f(x,y)','optimum iter')
grid on

----------------------------------------- RESPONDING TO THE COMMENT ----------------------------------------------------

clc
clear
close all

% defining anonymous functions for f,dx and dy for convenience in
% computation
f= @(x,y) -8.*x + 12.*y + x.^2 - 2.*x.^4 - 2.*x.*y + 4.*y.^2;
dx = @(x,y) -8+2*x-8*x^3-2*y;
dy = @(x,y) 12-2*x+8*y;

% S(j) indicates an array S,where j means position or index in the array
% S(1) refers to element at position 1 in array S
% S(j) = p means put value p at position j in the array S

% array can be taken to mean the same thing as vector

% assuming the vectors x and y both start at -4 --> at point (-4,-4).
% the starting point (-4,-4) is chosen because the function f seems to have an optimum point
% when both x and y are negative, starting at (0,0) would mean the optimisation may not
% converge. Try starting at (0,0) and view the response

x(1) = -4; y(1) = -4; %(-4,-4)

%definig the vector z to hold values of f
z(1) = f(x(1),y(1));

max_iter = 300; % chosen arbitrary, it can be any value greater than 0
for i=1:max_iter

x(i+1) = x(i) + 0.001*dx(x(i),y(i));
y(i+1) = y(i) + 0.001*dy(x(i),y(i));

z(i+1) = f(x(i+1),y(i+1));
diff(i) = z(i+1) - z(i);

end

optimum_iteration = find(diff==min(diff))
optimum_x = x(optimum_iteration)
optimum_y = y(optimum_iteration)
optimum_f = f(optimum_x,optimum_y)

% displaying results
figure
plot3(x,y,z)
hold on
plot3(optimum_x,optimum_y,optimum_f,'o g')
title('f = -8*x + 12*y + x^2 - 2*x^4 - 2*x*y + 4*y^2')
xlabel('x');
ylabel('y')
zlabel('f')
grid on
legend('f','optimum f')

figure
plot(diff,'r')
hold on
plot(optimum_iteration,diff(optimum_iteration),'x k')
title('Change in f with iteration')
ylabel('f_k_+_1 - f_k')
xlabel('iteration (k)')
legend('change in f','minimum change in f')
grid on

% combine diff and f vs iteration on one plot
% considerign less iterations for better viewing
iterations_considered = 1:100; % must be less than max_iter

figure
plot(z,'b')
grid on
title('f vs iterations')
xlabel('iteration (k)')
ylabel('f')

figure
plot(z(iterations_considered),'b')
hold on
plot(diff(iterations_considered),'r')
hold on
plot(optimum_iteration,optimum_f,'o g')
hold on
plot(optimum_iteration,diff(optimum_iteration),'x k')
title('Variation of f, change in f with iterations')
xlabel('iteration (k)')
legend('f','change in f','optimum f','minimum change in f')
grid on

--------------------------------------------- RESPONDING TO THE COMMENT----------------------------------

% initial/starting point for iterative methods like gradient ascent/descent is very important mostly when the objective %function is not so smooth (having several local optima)

YOU MIGHT FIND THE FOLLOWING HELPFUL:

clc
clear
close all

%using graphical approach to get a clue about the global maximum
% interval [-4,4] used for both x and y

f= @(x,y) -8.*x + 12.*y + x.^2 - 2.*x.^4 - 2.*x.*y + 4.*y.^2;
dx = @(x,y) -8+2*x-8*x^3-2*y;
dy = @(x,y) 12-2*x+8*y;

x_vector = linspace(-4,4,100); %[-4,4]
y_vector = linspace(-4,4,100); %[-4,4]
f_vector = f(x_vector,y_vector);
f_max = max(f_vector)
x_at_fmax = x_vector(find(f_vector==f_max,1))
y_at_fmax = y_vector(find(f_vector==f_max,1))

%{
figure
[X,Y] = meshgrid(x_vector,y_vector);
Z = f(X,Y);
mesh(X,Y,Z)
xlabel('x'),ylabel('y'),zlabel('f')
%}

figure
plot3(x_vector,y_vector,f_vector)
xlabel('x'),ylabel('y'),zlabel('f')
hold on
plot3(x_at_fmax,y_at_fmax,f_max,'x r')
grid on
title('Original plot of the f on interval [-4,4] for both x and y')

figure
plot(x_vector,f_vector)
xlabel('x'),ylabel('f')
hold on
plot(x_at_fmax,f_max,'x r')
grid on
title('2D view showing the actual global max value')

------------------------------------------------------------------------------------------------------------------------------------------

clc
clear
close all

f= @(x,y) -8.*x + 12.*y + x.^2 - 2.*x.^4 - 2.*x.*y + 4.*y.^2;
dx = @(x,y) -8+2*x-8*x^3-2*y;
dy = @(x,y) 12-2*x+8*y;

%starting at [-4,-4]
x0 = -4; y0 = -4;
x(1) = x0; y(1) = y0;
z(1) = f(x(1),y(1));

max_iter = 500;
for i=1:max_iter

x(i+1) = x(i) + 0.001*dx(x(i),y(i));
y(i+1) = y(i) + 0.001*dy(x(i),y(i));

z(i+1) = f(x(i+1),y(i+1));
end

plot(z)
title('f vs iteration')
xlabel('iteration (k)')
ylabel('f_k')
grid on

figure
plot(x,z)
title('f vs x')
xlabel('x')
ylabel('f')
grid on

figure
plot(y,z)
title('f vs y')
xlabel('y')
ylabel('f')
grid on
f_max = max(z)
x_at_fmax = x(find(z==f_max,1))
y_at_fmax = y(find(z==f_max,1))

%PLEASE NOTE THE OPTIMISATION HAS NOT CONVERGED TO THE VALUE EXPECTED
% ie f should converge to 5.1
  

----------------------------------------------------------------------------------------------------------------------------------------

clc
clear
close all

%let's start at [1,1] a point close to the known maximum value to minimiZe chances of getting stuck
%in local maxima

f= @(x,y) -8.*x + 12.*y + x.^2 - 2.*x.^4 - 2.*x.*y + 4.*y.^2;
dx = @(x,y) -8+2*x-8*x^3-2*y;
dy = @(x,y) 12-2*x+8*y;

x0 = 1; y0 = 1;
x(1) = x0; y(1) = y0;
z(1) = f(x(1),y(1));

max_iter = 500;
for i=1:max_iter

x(i+1) = x(i) + 0.001*dx(x(i),y(i));
y(i+1) = y(i) + 0.001*dy(x(i),y(i));

z(i+1) = f(x(i+1),y(i+1));
end

plot(z)
title('f vs iteration')
xlabel('iteration (k)')
ylabel('f_k')
grid on

figure
plot(x,z)
title('f vs x')
xlabel('x')
ylabel('f')
grid on

figure
plot(y,z)
title('f vs y')
xlabel('y')
ylabel('f')
grid on
f_max = max(z)
x_at_fmax = x(find(z==f_max,1))
y_at_fmax = y(find(z==f_max,1))

%PLEASE NOTE THE OPTIMISATION HAS NOT CONVERGED TO THE VALUE EXPECTED
% ie f should converge to 5.1
  

------------------------------------------------------------------------------------------------------------------------------

clc
clear
close all

% TWISTING THE GRADIENT ASCENT TO CONVERGE TO THE MAXIMUM VALUE

f= @(x,y) -8.*x + 12.*y + x.^2 - 2.*x.^4 - 2.*x.*y + 4.*y.^2;
dx = @(x,y) -8+2*x-8*x^3-2*y;
dy = @(x,y) 12-2*x+8*y;

% CONSIDER THE SIGNS OF GRADIENTS Dy, Dx at the initial iteration

%starting point
x0 = -4; y0 = -4;

Dx = sign(dx(x0,y0));
Dy = sign(dy(x0,y0));

% Since Dy is negative, yet positive gradient is needed to achieve maximum
% value, we need to force dy to the positive direction by using minus in
% the equation y(i+1) = y(i) - 0.001*dy(x(i),y(i))
% more like gradient descent but applied to one of the variables that tend
% to diverge from the positive course

x(1) = x0; y(1) = y0;
z(1) = f(x(1),y(1));

max_iter = 500;
for i=1:max_iter

x(i+1) = x(i) + 0.001*dx(x(i),y(i));
y(i+1) = y(i) - 0.001*dy(x(i),y(i));

z(i+1) = f(x(i+1),y(i+1));
end

f_max = max(z)
x_at_fmax = x(find(z==f_max,1))
y_at_fmax = y(find(z==f_max,1))

plot(z)
title('f vs iteration')
xlabel('iteration (k)')
ylabel('f_k')
grid on

figure
plot(x,z)
title('f vs x')
xlabel('x')
ylabel('f')
grid on
hold on
plot(x_at_fmax,f_max,'x r')

figure
plot(y,z)
title('f vs y')
xlabel('y')
ylabel('f')
grid on
hold on
plot(y_at_fmax,f_max,'x r')

figure
plot3(x,y,z)
hold on
plot3(x_at_fmax,y_at_fmax,f_max,'x r')
grid on

% f = 4.3485 is close enough to 5.1
% if a better starting point like [0.9,0.9] was chosen, we would be close to 5.1

COMMENT DOWN FOR ANY QUERY RELATED TO THIS ANSWER,
GIVE A THUMBS UP IF YOU LIKE THE ANSWER