From 44ca73fe54887fcdba4ef549c908b2d680f0936f Mon Sep 17 00:00:00 2001
From: mhjensen <morten.hjorth-jensen@fys.uio.no>
Date: Sat, 6 Oct 2018 13:26:00 +0200
Subject: [PATCH] added code for stochastic gradient descent

---
 doc/src/Splines/Splines.do.txt | 77 ++++++++++++++++++++++++++++++++++
 1 file changed, 77 insertions(+)

diff --git a/doc/src/Splines/Splines.do.txt b/doc/src/Splines/Splines.do.txt
index a4556fe71..ef15d68f9 100644
--- a/doc/src/Splines/Splines.do.txt
+++ b/doc/src/Splines/Splines.do.txt
@@ -1378,6 +1378,83 @@ print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
 !split
 =====  Program for stochastic gradient =====
 
+!bc pycod
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import SGDRegressor
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print("Own inversion")
+print(theta_linreg)
+sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
+sgdreg.fit(x,y.ravel())
+print("sgdreg from scikit")
+print(sgdreg.intercept_, sgdreg.coef_)
+
+
+theta = np.random.randn(2,1)
+
+eta = 0.1
+Niterations = 1000
+m = 100
+
+for iter in range(Niterations):
+    gradients = 2.0/m*xb.T.dot(xb.dot(theta)-y)
+    theta -= eta*gradients
+print("theta frm own gd")
+print(theta)
+
+xnew = np.array([[0],[2]])
+xbnew = np.c_[np.ones((2,1)), xnew]
+ypredict = xbnew.dot(theta)
+ypredict2 = xbnew.dot(theta_linreg)
+
+
+n_epochs = 50
+t0, t1 = 5, 50
+m = 100
+def learning_schedule(t):
+    return t0/(t+t1)
+
+theta = np.random.randn(2,1)
+
+for epoch in range(n_epochs):
+    for i in range(m):
+        random_index = np.random.randint(m)
+        xi = xb[random_index:random_index+1]
+        yi = y[random_index:random_index+1]
+        gradients = 2 * xi.T.dot(xi.dot(theta)-yi)
+        eta = learning_schedule(epoch*m+i)
+        theta = theta - eta*gradients
+print("theta from own sdg")
+print(theta)
+
+
+
+
+
+
+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()
+
+!ec
+
+
+
+
 !split
 ===== Momentum based methods =====