-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgradientvsnewton2D.m
127 lines (102 loc) · 3.29 KB
/
gradientvsnewton2D.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
clear all
close all
% define range
xvec = 0:5:180;
yvec = 0:5:180;
[X,Y] = meshgrid(xvec,yvec);
% define function
% NOTE: this is the misfit function (not the synthetic g(m))
fx = @(X,Y) (cosd(X).^2 + cosd(Y).^2).^2;
Z = fx(X,Y);
% plot function
figure(1)
surf(X,Y,Z);
title(sprintf('f(x,y) = (cos(x)^2 + cos(y)^2)^2'));
xlabel('x'); ylabel('y'); zlabel('f(x,y)')
%---------------------------------------------------------------
% calculate gradient
% ANALYTICALLY
% This can also be used in partial derivative matrix
dfdx = @(X,Y) -4*(cosd(X).*sind(X)).*(cosd(X).^2 + cosd(Y).^2);
dfdy = @(X,Y) -4*(cosd(Y).*sind(Y)).*(cosd(X).^2 + cosd(Y).^2);
% Plot gradient
figure
subplot(2,2,1)
quiver(X,Y,dfdx(X,Y),dfdy(X,Y),'k')
title('Analytical estimation of gradient');
xlabel('x'); ylabel('y');
%---------------------------------------------------------------
% Numerical estimation of gradient (without analytical derivatives)
delx = 5; dely = 5;
% remove boundary points
xvec_grad = xvec(2):delx:xvec(end-1);
yvec_grad = yvec(2):delx:yvec(end-1);
[numX,numY] = meshgrid(xvec_grad,yvec_grad);
% Numerical gradient using central difference
for ii=2:length(xvec)-1
for jj=2:length(yvec)-1
delzdelx(jj-1,ii-1) = (Z(jj,ii+1)-Z(jj,ii-1))/2*(delx);
delzdely(jj-1,ii-1) = (Z(jj+1,ii)-Z(jj-1,ii))/2*(dely);
end
end
subplot(2,2,2)
quiver(numX,numY,delzdelx,delzdely,'r')
title('Numerical estimation of gradient');
xlabel('x'); ylabel('y');
%---------------------------------------------------------------
% Compare with MATLAB inbuilt
[px,py] = gradient(Z);
subplot(2,2,3)
quiver(X,Y,px,py)
title('MATLAB in-built gradient');
xlabel('x'); ylabel('y');
%=================================================================
% Find minimum
% Gradient descent
% Hessian estimation not required
% Initial guess
x = 40;
y = 60;
delx = 1; dely = 1;
step_size = 10;
for ii=1:1000
% compute gradient analytically
dgdx = dfdx(x(ii),y(ii)); dgdy = dfdy(x(ii),y(ii));
% You can also compute gradient Numerically
%dgdx = (fx(x(ii)+delx,y(ii))- fx(x(ii)-delx,y(ii)))/(2*delx);
%dgdy = (fx(x(ii),y(ii)+dely)- fx(x(ii),y(ii)-dely))/(2*dely);
% Model update
x(ii+1) = x(ii) - step_size*dgdx;
y(ii+1) = y(ii) - step_size*dgdy;
end
figure(1)
hold on
plot3(x,y,fx(x,y),'o-r','LineWidth',2)
quiver(X,Y,px,py)
%=================================================================
% Find minimum
% Newtons method
% Hessian
d2fdx2 = @(X,Y) -4*(cosd(2*X)).*(cosd(X).^2 + cosd(Y).^2)+2*sind(2*X);
d2fdy2 = @(X,Y) -4*(cosd(2*Y)).*(cosd(X).^2 + cosd(Y).^2)+2*sind(2*Y);
d2fdxdy = @(X,Y) 2*(sind(2*X)).*sind(2*Y);
d2fdydx = @(X,Y) 2*(sind(2*X)).*sind(2*Y);
% Initial guess
x = 100;
y = 30;
for ii = 1:1500
Gh = [dfdx(x(ii),y(ii)) dfdy(x(ii),y(ii))];
% second order partial derivative matrix
K = [d2fdx2(x(ii),y(ii)) d2fdxdy(x(ii),y(ii));
d2fdydx(x(ii),y(ii)) d2fdy2(x(ii),y(ii))];
% Hessian gives the direction of descent
H = (Gh' * Gh) + K*fx(x(ii),y(ii));
dm = -inv(H) * Gh' * fx(x(ii),y(ii));
% Model update
x(ii+1) = x(ii) + dm(1);
y(ii+1) = y(ii) + dm(2);
end
figure(1)
hold on
plot3(x,y,fx(x,y),'o-k','LineWidth',2)
quiver(X,Y,px,py)