Morten Hjorth-Jensen
52 lines
1.1 KiB
Python
52 lines
1.1 KiB
Python
# 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()
|
|
|
|
|