diff --git a/doc/src/week40/codes/test1.py b/doc/src/week40/codes/test1.py new file mode 100644 index 000000000..a41ac37c0 --- /dev/null +++ b/doc/src/week40/codes/test1.py @@ -0,0 +1,51 @@ +# Using Autograd to calculate gradients +from random import random, seed +import numpy as np +import autograd.numpy as np +import matplotlib.pyplot as plt +from autograd import grad + +def CostOLS(beta): + return (1.0/n)*np.sum((y-X @ beta)**2) + +n = 100 +x = 2*np.random.rand(n,1) +y = 4+3*x+np.random.randn(n,1) + +X = np.c_[np.ones((n,1)), x] +XT_X = X.T @ X +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(theta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +theta = np.random.randn(2,1) +eta = 1.0/np.max(EigValues) +Niterations = 1000 + +training_gradient = grad(CostOLS) + +for iter in range(Niterations): + gradients = training_gradient(theta) + theta -= eta*gradients +print("theta from own gd") +print(theta) + +xnew = np.array([[0],[2]]) +Xnew = np.c_[np.ones((2,1)), xnew] +ypredict = Xnew.dot(theta) +ypredict2 = Xnew.dot(theta_linreg) + +plt.plot(xnew, ypredict, "r-") +plt.plot(xnew, ypredict2, "b-") +plt.plot(x, y ,'ro') +plt.axis([0,2.0,0, 15.0]) +plt.xlabel(r'$x$') +plt.ylabel(r'$y$') +plt.title(r'Random numbers ') +plt.show() + + diff --git a/doc/src/week40/codes/test2.py b/doc/src/week40/codes/test2.py new file mode 100644 index 000000000..2b46cc1b6 --- /dev/null +++ b/doc/src/week40/codes/test2.py @@ -0,0 +1,33 @@ +import autograd.numpy as np +from autograd import grad + +def sigmoid(x): + return 0.5 * (np.tanh(x / 2.) + 1) + +def logistic_predictions(weights, inputs): + # Outputs probability of a label being true according to logistic model. + return sigmoid(np.dot(inputs, weights)) + +def training_loss(weights): + # Training loss is the negative log-likelihood of the training labels. + preds = logistic_predictions(weights, inputs) + label_probabilities = preds * targets + (1 - preds) * (1 - targets) + return -np.sum(np.log(label_probabilities)) + +# Build a toy dataset. +inputs = np.array([[0.52, 1.12, 0.77], + [0.88, -1.08, 0.15], + [0.52, 0.06, -1.30], + [0.74, -2.49, 1.39]]) +targets = np.array([True, True, False, True]) + +# Define a function that returns gradients of training loss using Autograd. +training_gradient_fun = grad(training_loss) + +# Optimize weights using gradient descent. +weights = np.array([0.0, 0.0, 0.0]) +print("Initial loss:", training_loss(weights)) +for i in range(100): + weights -= training_gradient_fun(weights) * 0.01 + +print("Trained loss:", training_loss(weights)) diff --git a/doc/src/week40/codes/test3.py b/doc/src/week40/codes/test3.py new file mode 100644 index 000000000..306ff15f5 --- /dev/null +++ b/doc/src/week40/codes/test3.py @@ -0,0 +1,73 @@ +# Using Autograd to calculate gradients using SGD +# OLS example +from random import random, seed +import numpy as np +import autograd.numpy as np +import matplotlib.pyplot as plt +from autograd import grad + +def CostOLS(y,X,theta): + return np.sum((y-X @ theta)**2) + +n = 100 +x = 2*np.random.rand(n,1) +y = 4+3*x+np.random.randn(n,1) + +X = np.c_[np.ones((n,1)), x] +XT_X = X.T @ X +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(theta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +theta = np.random.randn(2,1) +eta = 1.0/np.max(EigValues) +Niterations = 1000 + +# Note that we request the derivative wrt third argument (theta, 2 here) +training_gradient = grad(CostOLS,2) + +for iter in range(Niterations): + gradients = (1.0/n)*training_gradient(y, X, theta) + theta -= eta*gradients +print("theta from own gd") +print(theta) + +xnew = np.array([[0],[2]]) +Xnew = np.c_[np.ones((2,1)), xnew] +ypredict = Xnew.dot(theta) +ypredict2 = Xnew.dot(theta_linreg) + +plt.plot(xnew, ypredict, "r-") +plt.plot(xnew, ypredict2, "b-") +plt.plot(x, y ,'ro') +plt.axis([0,2.0,0, 15.0]) +plt.xlabel(r'$x$') +plt.ylabel(r'$y$') +plt.title(r'Random numbers ') +plt.show() + +n_epochs = 50 +M = 5 #size of each minibatch +m = int(n/M) #number of minibatches +t0, t1 = 5, 50 +def learning_schedule(t): + return t0/(t+t1) + +theta = np.random.randn(2,1) + +for epoch in range(n_epochs): +# Can you figure out a better way of setting up the contributions to each batch? + for i in range(m): + random_index = np.random.randint(m) + xi = X[random_index*M:random_index*M+M] + yi = y[random_index*M:random_index*M+M] + gradients = (2.0/M)*training_gradient(yi, xi, theta) + eta = learning_schedule(epoch*m+i) + theta = theta - eta*gradients +print("theta from own sdg") +print(theta) +