Gradient Descent¶
In [1]:
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
Example¶
In [64]:
x, y = sp.symbols("x y", real = True)
f = x**4 + y**4 + sp.Rational(1, 2) * (x**2 + y**2) + 2 * x * y
f_num = sp.lambdify((x, y), f)
f_grad = sp.Matrix([sp.diff(f, x), sp.diff(f, y)])
f_grad_num = sp.lambdify((x, y), list(f_grad)) # use list here to get a nicer output as flat array
H = sp.hessian(f, (x, y))
H_num = sp.lambdify((x, y), H)
x_num = np.linspace(-1, 1, 400)
y_num = np.linspace(-2, 1, 400)
X_num, Y_num = np.meshgrid(x_num, y_num)
Z_num = f_num(X_num, Y_num)
plt.contour(X_num, Y_num, Z_num, 30)
Out[64]:
<matplotlib.contour.QuadContourSet at 0x783a7d77b5f0>
Exact Solutions¶
In [6]:
stationary_points = sp.solve([sp.diff(f, x), sp.diff(f, y)], [x, y], dict = True)
stationary_points
Out[6]:
[{x: -1/2, y: 1/2}, {x: 0, y: 0}, {x: 1/2, y: -1/2}]
In [16]:
eigenvalues_values = [H.subs({x: sol[x], y: sol[y]}).eigenvals() for sol in stationary_points]
eigenvalues_values
Out[16]:
[{6: 1, 2: 1}, {3: 1, -1: 1}, {6: 1, 2: 1}]
In [9]:
eigenvalues_values[0].keys() # eigenvalues are the keys of the list
Out[9]:
dict_keys([6, 2])
So we get the numerical value by ... and therefore saddle point for (-1/2, 1/2):
In [10]:
[ev.evalf() for ev in eigenvalues_values[0].keys()]
Out[10]:
[6.00000000000000, 2.00000000000000]
Saddlepoint for point for (0, 0):
In [11]:
[ev.evalf() for ev in eigenvalues_values[1].keys()]
Out[11]:
[3.00000000000000, -1.00000000000000]
Minimum for point for (1/2, -1/2):
In [13]:
[ev.evalf() for ev in eigenvalues_values[2].keys()]
Out[13]:
[6.00000000000000, 2.00000000000000]
Newton's Method¶
1 dimensional approach to solve g(x) = 0 (scalar equation)
sol_new = sol_old - g(sol_old) / g'(sol_old)
Approach to solve g(x, y) = [0, 0] (2 equations!)
sol_new, sol_old and g(sol_old) are vectors now
g'(sol_old) gets a matrix of all derivatives H(sol_old) = [[g1x, g1y], [g2x, g2y]](sol_old)
sol_new = sol_old - H^{-1}(sol_old) * g(sol_old)
Approach to solve grad f(x, y) = [0, 0] (2 equations!)
sol_new, sol_old and g = grad f(sol_old) are again vectors
g'(sol_old) gets a matrix of all derivatives H(sol_old) = [[g1x, g1y], [g2x, g2y]] = [[fxx, fxy], [fyx, fyy]](sol_old)
sol_new = sol_old - H^{-1}(sol_old) * g(sol_old)
In [39]:
sol_0 = np.array([0.99, 0.50])
g_0 = np.array(f_grad_num(sol_0[0], sol_0[1]))
H_0 = np.array(H_num(sol_0[0], sol_0[1]))
H_0_inv = np.linalg.inv(H_0)
sol_1 = sol_0 - H_0_inv @ g_0
sol_1
Out[39]:
array([ 0.61748733, -0.05874367])
In [40]:
g_1 = np.array(f_grad_num(sol_1[0], sol_1[1]))
H_1 = np.array(H_num(sol_1[0], sol_1[1]))
H_1_inv = np.linalg.inv(H_1)
sol_2 = sol_1 - H_1_inv @ g_1
sol_2
Out[40]:
array([ 1.08769541, -2.09044748])
In [41]:
g_2 = np.array(f_grad_num(sol_2[0], sol_2[1]))
H_2 = np.array(H_num(sol_2[0], sol_2[1]))
H_2_inv = np.linalg.inv(H_2)
sol_3 = sol_2 - H_2_inv @ g_2
sol_3
Out[41]:
array([ 0.86163572, -1.39980011])
In [66]:
sol_old = np.array([0.99, 0.50])
sol_new = np.array([0.99, 0.50])
solutions = [sol_old]
for i in range(50):
g_old = np.array(f_grad_num(sol_old[0], sol_old[1]))
H_old = np.array(H_num(sol_old[0], sol_old[1]))
# H_old_inv = np.linalg.inv(H_old)
# sol_new = sol_old - H_old_inv @ g_old
sol_new = sol_old - np.linalg.solve(H_old, g_old) # better to avoid inverse and ask numpy to solve a system of linear equations directly
solutions.append(sol_new)
sol_old = sol_new
sol_new
Out[66]:
array([ 0.5, -0.5])
In [67]:
path = np.array(solutions)
x_path = path[:, 0]
y_path = path[:, 1]
plt.contour(X_num, Y_num, Z_num, 30)
plt.plot(x_path, y_path, 'ro-', markersize=3)
Out[67]:
[<matplotlib.lines.Line2D at 0x783a7d689c70>]
We could improve the method, we see that the second solution is very bad, but then it is converging!
Gradient Descent¶
In [57]:
step_size = 0.1
sol_0 = np.array([0.99, 0.50])
sol_1 = sol_0 - step_size * np.array(f_grad_num(sol_0[0], sol_0[1]))
sol_1
Out[57]:
array([0.4028804, 0.202 ])
In [58]:
sol_2 = sol_1 - step_size * np.array(f_grad_num(sol_1[0], sol_1[1]))
sol_2
Out[58]:
array([0.29603533, 0.09792696])
In [59]:
sol_3 = sol_2 - step_size * np.array(f_grad_num(sol_2[0], sol_2[1]))
sol_3
Out[59]:
array([0.23646896, 0.02855156])
In [69]:
sol_old = np.array([0.99, 0.50])
sol_new = np.array([0.99, 0.50])
solutions = [sol_old]
for i in range(50):
sol_new = sol_old - step_size * np.array(f_grad_num(sol_old[0], sol_old[1]))
solutions.append(sol_new)
sol_old = sol_new
sol_new
Out[69]:
array([ 0.49970797, -0.49970797])
In [70]:
path = np.array(solutions)
x_path = path[:, 0]
y_path = path[:, 1]
plt.contour(X_num, Y_num, Z_num, 20)
plt.plot(x_path, y_path, 'ro-', markersize=3)
Out[70]:
[<matplotlib.lines.Line2D at 0x783a7d6c6960>]
