%matplotlib inline
import numpy as np
import scipy.optimize
import sklearn.datasets
import matplotlib.pyplot as plt
np.set_printoptions(suppress=True, precision=6, linewidth=200)
plt.style.use('ggplot')
def f(x, y):
return x ** 2 + y ** 2 + x * (y + 2) + np.cos(3 * x)
def grad_x_f(x, y):
return 2 * x - 3 * np.sin(3 * x) + y + 2
def grad_y_f(x, y):
return x + 2 * y
def plot_f_contours():
xx, yy = np.meshgrid(np.linspace(-5, 5), np.linspace(-5, 5))
zz = f(xx, yy)
plt.contourf(xx, yy, zz, 50)
plt.contour(xx, yy, zz, 50, alpha=0.2, colors='black', linestyles='solid')
plt.xlabel('x')
plt.ylabel('y')
plt.figure(figsize=(10, 7))
plot_f_contours()
def optimize_f(x, y, step_size, steps):
# keep track of the parameters we tried so far
x_hist, y_hist = [x], [y]
# run gradient descent for the number of steps
for step in range(steps):
# compute the gradients at the current point
dx = grad_x_f(x, y)
dy = grad_y_f(x, y)
# apply the gradient descent updates to x and y
x = x - (step_size * dx)
y = y - (step_size * dy)
# store the new parameters
x_hist.append(x)
y_hist.append(y)
return x, y, f(x, y), x_hist, y_hist
# helper function that plots the results of the gradient descent optimization
def plot_gradient_descent_results(x, y, val, x_hist, y_hist):
# plot the path on the contour plot
plt.figure(figsize=(20, 7))
plt.subplot(1, 2, 1)
plot_f_contours()
plt.plot(x_hist, y_hist, '.-')
# plot the learning curve
plt.subplot(1, 2, 2)
plt.plot(f(np.array(x_hist), np.array(y_hist)), '.r-')
plt.title('Minimum value: %f' % f(x_hist[-1], y_hist[-1]))
results = optimize_f(x=3, y=2, step_size=0.1, steps=10)
plot_gradient_descent_results(*results)
# TODO: tune the parameters to find a better optimum
results = optimize_f(x=3, y=2, step_size=0.16, steps=100)
plot_gradient_descent_results(*results)
def optimize_f(x, y, step_size, steps, decay=1.0):
# keep track of the parameters we tried so far
x_hist, y_hist = [x], [y]
# run gradient descent for the number of steps
for step in range(steps):
# compute the gradients at this point
dx = grad_x_f(x, y)
dy = grad_y_f(x, y)
# apply the gradient descent updates to x and y
x = x - (step_size * ( decay ** step ) * dx) # TODO: compute the update including step size decay
y = y - (step_size * ( decay ** step ) * dy) # TODO: compute the update including step size decay
# store the new parameters
x_hist.append(x)
y_hist.append(y)
return x, y, f(x, y), x_hist, y_hist
# TODO: tune the parameters to find the local optimum
results = optimize_f(x=3, y=2, step_size=0.16, steps=100)
plot_gradient_descent_results(*results)
def sigmoid(x):
sig = 1 / (1 + np.exp(-x)) # TODO: implement the sigmoid function
return sig
def sigmoid_grad(x):
grad = (sigmoid(x) * (1 - sigmoid(x))) # TODO: implement the gradient of the sigmoid function
return grad
# try with a random input
x = np.random.uniform(-10, 10, size=5)
print('x:', x)
print('sigmoid(x):', sigmoid(x))
print('sigmoid_grad(x):', sigmoid_grad(x))
# start with some random inputs
x = np.random.uniform(-2, 2, size=5)
# compute the symbolic gradient
print('Symbolic', sigmoid_grad(x))
# TODO: compute the numerical gradient
def sig_numerical(x):
e = 0.0001
num = (sigmoid(x+(0.5*e)) - sigmoid(x-(0.5*e)))/e
return num
print('Numerical', sig_numerical(x))
def relu(x):
return np.maximum(x,0)
def relu_grad(x):
return np.greater(x,0).astype(int)
# try with a random input
x = np.random.uniform(-10, 10, size=5)
print('x:', x)
print('relu(x):', relu(x))
print('relu_grad(x):', relu_grad(x))
print()
# TODO: compute and compare the symbolic and numerical gradients
def relu_numerical(x):
e = 0.0001
num = (relu(x+(0.5*e)) - relu(x-(0.5*e)))/e
return num
print('Numerical', relu_numerical(x))
x = np.linspace(-10, 10, 100)
plt.figure(figsize=(15, 8))
plt.subplot(2, 2, 1)
plt.plot(x, sigmoid(x), label='Sigmoid')
plt.xlabel('x')
plt.legend(loc='upper left')
plt.subplot(2, 2, 2)
plt.plot(x, relu(x), label='ReLU')
plt.xlabel('x')
plt.legend(loc='upper left')
plt.subplot(2, 2, 3)
plt.plot(x, sigmoid_grad(x), label='Sigmoid gradient')
plt.xlabel('x')
plt.legend(loc='upper left')
plt.subplot(2, 2, 4)
plt.plot(x, relu_grad(x), label='ReLU gradient')
plt.xlabel('x')
plt.legend(loc='upper left');
def bce_loss(y, y_hat):
#print('vals:', y)
#print('predictions:', y_hat)
return -( y * np.log(y_hat) + (1-y) * np.log(1-y_hat) )
def bce_loss_grad(y, y_hat):
return (y_hat - y) / (y_hat - y_hat ** 2)
# try with some random inputs
y = np.random.randint(2, size=5)
y_hat = np.random.uniform(0, 1, size=5)
print('y:', y)
print('y_hat:', y_hat)
print('bceloss(y, y_hat):', bce_loss(y, y_hat))
print('bceloss_grad(y, y_hat):', bce_loss_grad(y, y_hat))
# TODO: compute and compare the symbolic and numerical gradients
def numerical_gradient_bce(y, y_hat):
e = 0.0001
y_hat_increment = y_hat + e
y_hat_decrement = y_hat - e
return ( bce_loss(y, y_hat_increment) - bce_loss(y, y_hat_decrement) ) / (2 * e)
print('numerical gradient:', numerical_gradient_bce(y, y_hat))
# initialize parameters
w = np.random.uniform(size=5)
b = np.random.rand()
# implement the model
def fn(x, y):
#print(x.shape)
# TODO: forward: compute h, y_hat, loss
h = np.matmul(np.transpose(x), w) + b
#print('h:', h)
y_hat = sigmoid(h)
loss = bce_loss(y, y_hat)
# TODO: backward: compute grad_y_hat, grad_h, grad_x
grad_y_hat = bce_loss_grad(y, y_hat)
#print(grad_y_hat)
grad_h = grad_y_hat * sigmoid_grad(h)
grad_x = grad_h * w
return loss, grad_x
# test with a random input
x = np.random.uniform(size=5)
y = 1
loss, grad_x = fn(x, y)
print("Loss", loss)
print("Gradient", grad_x)
# start with some random inputs
x = np.random.uniform(size=5)
y = 1
# set epsilon to a small value
eps = 0.00001
numerical_grad = np.zeros(x.shape)
# compute the gradient for each element of x separately
for i in range(len(x)):
# compute inputs at -eps/2 and +eps/2
x_a, x_b = x.copy(), x.copy()
x_a[i] += eps / 2
x_b[i] -= eps / 2
# compute the gradient for this element
loss_a, _ = fn(x_a, y)
loss_b, _ = fn(x_b, y)
numerical_grad[i] = (loss_a - loss_b) / eps
# compute the symbolic gradient
loss, symbolic_grad = fn(x, y)
print("Symbolic gradient")
print(symbolic_grad)
print("Numerical gradient")
print(numerical_grad)
# Computes y = x * w + b.
class Linear:
def __init__(self, n_in, n_out):
# initialize the weights randomly,
# using the Xavier initialization rule for scale
a = np.sqrt(6 / (n_in * n_out))
self.W = np.random.uniform(-a, a, size=(n_in, n_out))
self.b = np.zeros((n_out,))
def forward(self, x):
# TODO: compute the forward pass
y = np.matmul(x, self.W) + self.b
return y
def backward(self, x, dy):
# TODO: compute the backward pass,
# given dy, compute the gradients for x, W and b
dx = np.matmul(dy, self.W.transpose())
#print('x:', x.shape, 'dx:', dx.shape)
self.dW = np.matmul(x.transpose(), dy)
#print('W:', self.W)
self.db = np.sum(dy, axis=0)
#print('b:', self.b)
return dx
def step(self, step_size):
self.W = self.W - step_size * self.dW
self.b = self.b - step_size * self.db
def __str__(self):
return 'Linear %dx%d' % self.W.shape
# Try the new class with some random values.
# Debugging tip: always choose a unique length for each
# dimension, so you'll get an error if you mix them up.
x = np.random.uniform(size=(3, 5))
layer = Linear(5, 7)
y = layer.forward(x)
dx = layer.backward(x, np.ones_like(y))
print('x:', x)
print('y:', y)
print('dx:', dx)
# Computes y = 1 / (1 + exp(-x)).
class Sigmoid:
def forward(self, x):
return sigmoid(x)
def backward(self, x, dy):
# TODO: compute the backward pass,
return dy * sigmoid_grad(x)
def step(self, step_size):
pass # ? is something supposed to happen here
def __str__(self):
return 'Sigmoid'
# try the new class with some random values
x = np.random.uniform(size=(3, 5))
layer = Sigmoid()
y = layer.forward(x)
dx = layer.backward(x, np.ones_like(y))
print('y:', y)
print('dx:', dx)
# Computes y = max(0, x).
class ReLU:
def forward(self, x):
return relu(x)
def backward(self, x, dy):
# TODO: compute the backward pass,
return dy * relu_grad(x)
def step(self, step_size):
pass # ? Might add something here later
def __str__(self):
return 'ReLU'
# try the new class with some random values
x = np.random.uniform(-10, 10, size=(3, 5))
layer = ReLU()
y = layer.forward(x)
dx = layer.backward(x, np.ones_like(y))
print('y:', y)
print('dx:', dx)
## Verify gradient computations for Linear
# test for dx
layer = Linear(5, 7)
def test_fn(x):
x = x.reshape(3, 5)
# multiply the output with a constant to check if
# the gradient uses dy
return 2 * np.sum(layer.forward(x))
def test_fn_grad(x):
x = x.reshape(3, 5)
# multiply the incoming dy gradient with a constant
return layer.backward(x, 2 * np.ones((3, 7))).flatten()
err = scipy.optimize.check_grad(test_fn, test_fn_grad,
np.random.uniform(-10, 10, size=3 * 5))
print("err on dx:", "OK" if np.abs(err) < 1e-5 else "ERROR", err)
# test for dW
x = np.random.uniform(size=(3, 5))
layer = Linear(5, 7)
def test_fn(w):
layer.W = w.reshape(5, 7)
# multiply the output with a constant to check if
# the gradient uses dy
return 2 * np.sum(layer.forward(x))
def test_fn_grad(w):
layer.W = w.reshape(5, 7)
# multiply the incoming dy gradient with a constant
layer.backward(x, 2 * np.ones((3, 7)))
return layer.dW.flatten()
err = scipy.optimize.check_grad(test_fn, test_fn_grad,
np.random.uniform(-10, 10, size=5 * 7))
print("err on dW:", "OK" if np.abs(err) < 1e-5 else "ERROR", err)
# test for db
x = np.random.uniform(size=(3, 5,))
layer = Linear(5, 7)
def test_fn(b):
layer.b = b
# multiply the output with a constant to check if
# the gradient uses dy
return 2 * np.sum(layer.forward(x))
def test_fn_grad(b):
layer.b = b
# multiply the incoming dy gradient with a constant
layer.backward(x, 2 * np.ones((x.shape[0], 7)))
return layer.db
err = scipy.optimize.check_grad(test_fn, test_fn_grad,
np.random.uniform(-10, 10, size=7))
print("err on db:", "OK" if np.abs(err) < 1e-5 else "ERROR", err)
## Verify gradient computation for Sigmoid
# test for dx
layer = Sigmoid()
def test_fn(x):
# multiply the output with a constant to check if
# the gradient uses dy
return np.sum(2 * layer.forward(x))
def test_fn_grad(x):
# multiply the incoming dy gradient with a constant
return layer.backward(x, 2 * np.ones(x.shape))
err = scipy.optimize.check_grad(test_fn, test_fn_grad,
np.random.uniform(-10, 10, size=5))
print("err on dx:", "OK" if np.abs(err) < 1e-5 else "ERROR", err)
## Verify gradient computation for ReLU
# test for dx
layer = ReLU()
def test_fn(x):
# multiply the output with a constant to check if
# the gradient uses dy
return 2 * np.sum(layer.forward(x))
def test_fn_grad(x):
# multiply the incoming dy gradient with a constant
return layer.backward(x, 2 * np.ones(x.shape))
err = scipy.optimize.check_grad(test_fn, test_fn_grad,
np.random.uniform(1, 10, size=5))
print("err on dx:", "OK" if np.abs(err) < 1e-5 else "ERROR", err)
class Net:
def __init__(self, layers):
self.layers = layers
def forward(self, x):
# compute the forward pass for each layer
trace = []
for layer in self.layers:
# compute the forward pass
y = layer.forward(x)
# store the original input for the backward pass
trace.append((layer, x))
x = y
# return the final output and the history trace
return y, trace
def backward(self, trace, dy):
# compute the backward pass for each layer
for layer, x in trace[::-1]:
# compute the backward pass using the original input x
dy = layer.backward(x, dy)
def step(self, learning_rate):
# apply the gradient descent updates of each layer
for layer in self.layers:
layer.step(learning_rate)
def __str__(self):
return '\n'.join(str(l) for l in self.layers)
# load the first two classes of the digits dataset
dataset = sklearn.datasets.load_digits()
digits_x, digits_y = dataset['data'], dataset['target']
# create a binary classification problem
digits_y = (digits_y < 5).astype(float)
# plot some of the digits
plt.figure(figsize=(10, 2))
plt.imshow(np.hstack([digits_x[i].reshape(8, 8) for i in range(10)]), cmap='gray')
plt.grid(False)
plt.tight_layout()
plt.axis('off')
# normalize the values to [0, 1]
digits_x -= np.mean(digits_x)
digits_x /= np.std(digits_x)
# print some statistics
print('digits_x.shape:', digits_x.shape)
print('digits_y.shape:', digits_y.shape)
print('min, max values:', np.min(digits_x), np.max(digits_x))
print('labels:', np.unique(digits_y))
# make a 50%/50% train/test split
train_prop = 0.5
n_train = int(digits_x.shape[0] * train_prop)
# shuffle the images
idxs = np.random.permutation(digits_x.shape[0])
# take a subset
x = {'train': digits_x[idxs[:n_train]],
'test': digits_x[idxs[n_train:]]}
y = {'train': digits_y[idxs[:n_train]],
'test': digits_y[idxs[n_train:]]}
print('Training samples:', x['train'].shape[0])
print('Test samples:', x['test'].shape[0])
def fit(net, x, y, epochs=25, learning_rate=0.001, mb_size=10):
# initialize the loss and accuracy history
loss_hist = {'train': [], 'test': []}
accuracy_hist = {'train': [], 'test': []}
for epoch in range(epochs):
# initialize the loss and accuracy for this epoch
loss = {'train': 0.0, 'test': 0.0}
accuracy = {'train': 0.0, 'test': 0.0}
# first train on training data, then evaluate on the test data
for phase in ('train', 'test'):
# compute the number of minibatches
steps = x[phase].shape[0] // mb_size
# loop over all minibatches
for step in range(steps):
# get the samples for the current minibatch
x_mb = x[phase][(step * mb_size):((step + 1) * mb_size)]
y_mb = y[phase][(step * mb_size):((step + 1) * mb_size), None]
# compute the forward pass through the network
pred_y, trace = net.forward(x_mb)
# compute the current loss and accuracy
loss[phase] += np.mean(bce_loss(y_mb, pred_y))
accuracy[phase] += np.mean((y_mb > 0.5) == (pred_y > 0.5))
# only update the network in the training phase
if phase == 'train':
# compute the gradient for the loss
dy = bce_loss_grad(y_mb, pred_y)
# backpropagate the gradient through the network
net.backward(trace, dy)
# update the weights
net.step(learning_rate)
# compute the mean loss and accuracy over all minibatches
loss[phase] = loss[phase] / steps
accuracy[phase] = accuracy[phase] / steps
# add statistics to history
loss_hist[phase].append(loss[phase])
accuracy_hist[phase].append(accuracy[phase])
print('Epoch %3d: loss[train]=%7.4f accuracy[train]=%7.4f loss[test]=%7.4f accuracy[test]=%7.4f' %
(epoch, loss['train'], accuracy['train'], loss['test'], accuracy['test']))
# plot the learning curves
plt.figure(figsize=(20, 5))
plt.subplot(1, 2, 1)
for phase in loss_hist:
plt.plot(loss_hist[phase], label=phase)
plt.title('BCE loss')
plt.xlabel('Epoch')
plt.legend()
plt.subplot(1, 2, 2)
for phase in accuracy_hist:
plt.plot(accuracy_hist[phase], label=phase)
plt.title('Accuracy')
plt.xlabel('Epoch')
plt.legend()
# construct network
net = Net([
Linear(64, 32),
ReLU(),
Linear(32, 1),
Sigmoid()])
# TODO: tune the hyperparameters
fit(net, x, y,
epochs = 25,
learning_rate = 0.01,
mb_size = 32)
# construct network
net = Net([
Linear(64, 32),
Linear(32, 1),
Sigmoid()])
# TODO: tune the hyperparameters
fit(net, x, y,
epochs = 25,
learning_rate = 0.01,
mb_size = 32)
# construct network
net = Net([
Linear(64, 1),
Sigmoid()])
# TODO: tune the hyperparameters
fit(net, x, y,
epochs = 25,
learning_rate = 0.01,
mb_size = 32)
net = Net([
Linear(64, 32),
ReLU(),
Linear(32, 16),
ReLU(),
Linear(16, 8),
ReLU(),
Linear(8, 4),
ReLU(),
Linear(4, 1),
Sigmoid()])
# TODO: tune the hyperparameters
fit(net, x, y,
epochs = 25,
learning_rate = 0.01,
mb_size = 32)