<jupyter_start><jupyter_text>Distributed optimization algorithmsThe goal of this notebook is to work on distributed optimization algorithms, which are the foundation for large scale analytics and machine learning. Specifically, we will focus on the details of stochastic gradient descent (SGD). To do so, we will work on a simple regression problem, where we will apply SGD to minimize a loss function, as defined for the problem at hand. The emphasis of this laboratory is **not** on the machine learning part: even if you've never worked on regression problems, this shouldn't prevent you from being successful in developing the Notebook.Next, an outline of the steps we will follow in this Notebook:* Breif introduction to linear regression* Implementation of serial algorithms: from Gradient Descent, to Stochastic Gradient Descent* Implementation of distributed algorithms with Apache Spark  Initialization code<jupyter_code>%matplotlib inline
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy import stats 

from sklearn.datasets.samples_generator import make_regression<jupyter_output>/opt/conda/lib/python3.5/site-packages/sklearn/utils/fixes.py:64: DeprecationWarning: inspect.getargspec() is deprecated, use inspect.signature() instead
  if 'order' in inspect.getargspec(np.copy)[0]:<jupyter_text>A simple example: linear regressionLet’s see briefly how to use gradient descent in a simple least squares regression setting. Asssume we have an output variable $y$ which we think depends linearly on the input vector $x$. That is, we have:$$x = \begin{bmatrix}x_1 \\ x_2 \\ \cdots \\ x_m\end{bmatrix}, y = \begin{bmatrix}y_1 \\ y_2 \\ \cdots \\ y_m\end{bmatrix}$$We approximate $y_i$ by:$$ f_{\theta}(x_i) = \theta_1 + \theta_2 x_i$$Define the loss function for the simple linear least squares regression as follows:$$ J(\theta) = \frac{1}{2} \sum_{i=1}^{m}(f_{\theta}(x_i) - y_i)^2$$ Now, let's use scikit learn to create a regression problem. A few notes are in order:* The call to ```make_regression``` essentially generates samples for a regression problem* The call to ```stats.linegress``` calculates a linear least-squares regression for two sets of measurementsThis means we have a sort of "basline" to experiment with our SGD implementation.<jupyter_code>x, y = make_regression(n_samples = 10000, 
                       n_features=1, 
                       n_informative=1, 
                       noise=20,
                       random_state=2017)
x = x.flatten()
slope, intercept,_,_,_ = stats.linregress(x,y)
best_fit = np.vectorize(lambda x: x * slope + intercept)
plt.plot(x,y, 'o', alpha=0.5)
grid = np.arange(-3,3,0.1)
plt.plot(grid,best_fit(grid), '.')<jupyter_output><empty_output><jupyter_text>Batch gradient descentBefore delving into SGD, let's take a simpler approach. Assume that we have a vector of paramters $\theta$ and a loss function $J(\theta)$, which we want to minimize. The loss function we defined above has the form:$$ J(\theta) = \sum_{i=1}^{m} J_i(\theta)$$where $J_i$ is associated with the i-th observation in our data set, such as the one we generated above. The batch gradient descent algorithm, starts with some initial feasible parameter $\theta$ (which we can either fix or assign randomly) and then repeatedly performs the update:$$ \theta^{(t+1)} = \theta^{(t)} - \eta \nabla_{\theta} J(\theta^{(t)}) = \theta^{(t)} - \eta \sum_{i=1}^{m} \nabla_{\theta} J_i(\theta^{(t)}) $$where $t$ is an iteration index, and $\eta$ is a constant controlling step-size and is called the learning rate. Note that in order to make a **single update**, we need to calculate the gradient **using the entire dataset**. This can be very inefficient for large datasets, and it is the goal of this Notebook to insist on this aspect.In code, the main loop for batch gradient descent looks like this:```pythonfor i in range(n_epochs):  params_grad = evaluate_gradient(loss_function, data, params)  params = params - learning_rate * params_grad```For a given number of iterations (also called epochs) $n_{e}$ , we first evaluate the gradient vector of the loss function using **ALL** examples in the data set, and then we update the parameters with a given learning rate. Batch gradient descent is guaranteed to converge to the global minimum for convex loss surfaces and to a local minimum for non-convex surfaces. NOTE: Who computes the gradient?Given a loss function $J(\theta)$, the gradient with respect to paramters $\theta$ must be derived manually. In other words, given the expression of $J(\theta)$, pencil and paper are required to derive the analytical form of its gradient. Then this expression can be plugged into our code.Recently, machine learning libraries have adopted the techniques of **automatic differentiation**, which eliminate this tedious and error prone step. Given a loss function $J(\theta)$, such libraries automatically compute the gradient. For the linear regression case, let's derive the update step for gradient descent. Recall that we have defined:$$ J(\theta) = \frac{1}{2} \sum_{i=1}^{m}(f_{\theta}(x_i) - y_i)^2$$$$f_\theta(x_i) = \theta_1 + \theta_2 x_i$$So we have that:$$ \nabla_\theta J(\theta) = \frac{1}{2} \sum_{i=1}^{m} \nabla_\theta (\theta_1 + \theta_2 x_i - y_i)^2 $$If we explicit the partial derivatives of the gradient, we have:$$\frac{\partial J(\theta)}{\partial \theta_1} = \sum_{i=1}^{m} (\theta_1 + \theta_2 x_i - y_i)$$$$\frac{\partial J(\theta)}{\partial \theta_2} = \sum_{i=1}^{m} (\theta_1 + \theta_2 x_i - y_i) x_i$$So now we can explicit the update rules for the two model parameters:$$ \theta_1^{(t+1)} = \theta_1^{(t)} - \eta \sum_{i=1}^{m}  (\theta_1^{(t)} + \theta_2^{(t)} x_i - y_i)$$$$ \theta_2^{(t+1)} = \theta_2^{(t)} - \eta \sum_{i=1}^{m}  (\theta_1^{(t)} + \theta_2^{(t)} x_i - y_i) x_i$$  Using matrix notationNow, before writing some code, let's see how can we simplify the above expressions using matrices. Note that this is not only useful for working on a more compact notation, but it helps reason about efficient computations using libraries such as numpy, which we will use extensively.We firstly expressed our prediction as: $f_\theta (x_i) = \theta_1 + \theta_2 x_i$. Let's introduce a surrogate dimension for our input set $x$, such that:$$x = \begin{bmatrix} 1& x_{1} \\ 1 & x_{2} \\ \cdots & \cdots \\ 1 & x_{m}\end{bmatrix}$$where we define, with a small abuse of notation, $x_i = \begin{bmatrix} 1 & x_{i}\end{bmatrix}$. Also, let's define vector $\theta = \begin{bmatrix} \theta_1 \\ \theta_2\end{bmatrix}$. Then, we can rewrite $f_\theta(x_i) = x_i \theta = \begin{bmatrix} 1 & x_{i}\end{bmatrix} \begin{bmatrix} \theta_1 \\ \theta_2\end{bmatrix} = \theta_1 + \theta_2 x_i$.Let's use this notation to rewrite our gradients in matrix form.$\nabla_\theta J(\theta) =  \frac{1}{2} \sum_{i=1}^{m} \nabla_\theta (x_i \theta - y_i)^2$$= \sum_{i=1}^{m} (x_i \theta - y_i) x_i^T =$$= x^T (x \theta - y) =$$= \begin{bmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_m \end{bmatrix} \Bigg( \begin{bmatrix} 1& x_{1} \\ 1 & x_{2} \\ \cdots & \cdots \\ 1 & x_{m}\end{bmatrix} \begin{bmatrix} \theta_1 \\ \theta_2\end{bmatrix} - \begin{bmatrix} y_1 \\ y_2 \\ \cdots \\ y_m\end{bmatrix} \Bigg)$  Numpy arraysWith the work we did above, we can now cast everything into numpy arrays, which are efficient, and for which an efficient implementation of vector and matrix operations exists. Specifically, above we used the traditional matrix notation, where we manipulate column vectors. Hence, we express matrix operations (namely matrix products) using the traditional "row-by-column" approach.In numpy, we avoid this formalism by using dot product operations. So, given two column vectors:$$a= \begin{bmatrix} a_1 \\ a_2 \\ \cdots \\ a_{m}\end{bmatrix} b= \begin{bmatrix} b_1 \\ b_2 \\ \cdots \\ b_{m}\end{bmatrix}$$we have that $a b^T = a \cdot b = a_1 b_1 + a_2 b_2 + \cdots + a_m b_m$, where $\cdot$ is the symbol we use for dot product. Question 1. Implement your own version of Gradient Descent, as a serial algorithm.Follow the guidelines below:     Define a function to perform gradient descent. The function should accept as inputs: the training data $x$ and $y$, the initial guess for the parameters ($\theta_1$ and $\theta_2$), the learning rate. Additional arguments include the definition of the maximum number of iterations before the algorithm stops, and a second stop condition on the marginal improvment on the loss.     Keep track of the values of the loss, for each iteration.     Keep track of the gradient values, for each iteration.Once the ```gradient_descent``` function is defined, you can generate input data according to the cell above, that use scikitlearn.The output of your cell should contain the following information:     The values of the paramters obtained through Gradient Descent optimization      The values of the paramters obtained with the above cell, using scikitlearn      A plot of the loss versus iterations      A plot of the path the gradient takes from its initial to its final position<jupyter_code>def reshape(x):
    return np.c_[np.ones((x.shape[0])), x]


def loss(x,y,theta):
    lossRes=0
    for i in range (len(x)):
        lossRes+=(theta[0] + theta[1]*x[i]-y[i])**2
    return(lossRes/(2*len(x)))
    
def gradient(x,y,theta):
    res=np.dot(x.transpose(),(np.dot(x,theta)-y))
    return(res/len(x))

def norm(t1,t2):
    return((t1[0]-t2[0])**2+(t1[1]-t2[1])**2)



def gradient_descent(x,Y,theta1, theta2,learning_rate, n_iteration, marginal_loss):
    
    X = reshape(x)
    theta=np.array([theta1, theta2])
    
    gradientRes = gradient(X,Y,theta)
    theta=theta.transpose()
    thetaNew=theta - learning_rate*gradientRes
    print("theta",theta)
    
    lossRes=loss(x,Y,theta)
    
    lossHistory=[lossRes]
    gradientHistory=[]
    
    T1=[theta1]
    T2=[theta2]
    
    iteration=0
    

    
    while (iteration< n_iteration and norm(theta,thetaNew)>marginal_loss) : 
        
        gradientRes = gradient(X,Y,theta)
        gradientHistory.append(gradientRes)
        
        
        theta, thetaNew= thetaNew, thetaNew - learning_rate*gradientRes
        print("theta","thetaNew",theta,thetaNew)
        
        T1.append(theta[0])
        T2.append(theta[1])
        
        print("loss",loss(x,Y,theta))
        lossHistory.append(loss(x,Y,thetaNew))
        
        iteration+=1
        
        
    return(lossHistory,gradientHistory,theta,T1,T2)

n_iteration=100
(lossHistory,gradientHistory,theta,T1,T2) = gradient_descent(x,y,1, 1,0.05, n_iteration, 0.0001)

print(theta)
print(type(lossHistory),type(gradientHistory))
print(len(lossHistory),len(gradientHistory))

abscisses=[i for i in range (len(lossHistory))]
print(type(abscisses))

plt.subplot(2, 1, 1)
plt.plot(abscisses, lossHistory)


print("abs: {} gradh:{}".format(abscisses, gradientHistory))
gradientX=[gradientHistory[i][0] for i in range(len(gradientHistory))]
gradientY=[gradientHistory[i][1] for i in range(len(gradientHistory))]
plt.subplot(2, 1, 2)

plt.plot(gradientX, gradientY)

plt.show()<jupyter_output>theta [1 1]
theta thetaNew [ 0.93684426  1.66009444] [ 0.87368851  2.32018889]
loss 282.841790076
theta thetaNew [ 0.87368851  2.32018889] [ 0.81407223  2.94803839]
loss 274.692770843
theta thetaNew [ 0.81407223  2.94803839] [ 0.7579954   3.54364296]
loss 267.34146029
theta thetaNew [ 0.7579954   3.54364296] [ 0.70526241  4.108578  ]
loss 260.727812053
theta thetaNew [ 0.70526241  4.108578  ] [ 0.65567765  4.64441888]
loss 254.778671854
theta thetaNew [ 0.65567765  4.64441888] [ 0.60905619  5.15266403]
loss 249.427302498
theta thetaNew [ 0.60905619  5.15266403] [ 0.56522379  5.63473488]
loss 244.613629662
theta thetaNew [ 0.56522379  5.63473488] [ 0.52401634  6.09197964]
loss 240.283615582
theta thetaNew [ 0.52401634  6.09197964] [ 0.48527924  6.52567705]
loss 236.388654246
theta thetaNew [ 0.48527924  6.52567705] [ 0.4488669  6.93704  ]
loss 232.885025388
theta thetaNew [ 0.4488669  6.93704  ] [ 0.41464218  7.32721886]
loss 229.733403283
theta thetaNew [ 0.41464218  7.32721886] [ 0.38[...]<jupyter_text>Question 2. Plot the regression line, along with the training data, given the coefficients $\theta$ that you have obtained with Gradient Descent.<jupyter_code>X= np.c_[np.ones((x.shape[0])), x]
Y=X.dot(theta)

plt.figure()
plt.scatter(x, y, marker="o", c=y)
plt.plot(x, Y,"r-")<jupyter_output><empty_output><jupyter_text>Question 3. Plot a 3D surface representing: on the x,y axes the parameter values, on the z axis the loss value. Additionally, plot the trajectory of the loss function on the 3D surface, using the history you collected in the gradient_descent function you designed.Finally, plot a contour projection of the 3D surface, along with the corresponding projection of the trajectory followed by your Gradient Descent algorithm.<jupyter_code>from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

print(len(T1),len(T2),len(lossHistory))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(T1,T2,lossHistory)
ax.plot(T1,T2,lossHistory)
plt.show()
#map=cm.coolwarm,linewidth=0, antialiased=False)<jupyter_output>82 82 82<jupyter_text>Question 4. Plot the path the gradient takes from its initial to its final position.This is a two dimensional plot (because our parameter vector has size 2), with a point for each gradient value, and a line connecting the points.<jupyter_code>gradientX=[gradientHistory[i][0] for i in range(len(gradientHistory))]
gradientY=[gradientHistory[i][1] for i in range(len(gradientHistory))]

plt.plot(gradientX, gradientY)
plt.ylabel('gradientY')
plt.xlabel('gradientX')

plt.show()<jupyter_output><empty_output><jupyter_text>Stochastic Gradient DescentThe gradient descent algorithm makes intuitive sense as it always proceeds in the direction of steepest descent (the gradient of $J$) and guarantees that we find a local minimum (global under certain assumptions on $J$). When we have very large data sets, however, computing $\nabla_\theta J(\theta)$ can be computationally challenging: as noted above, we must process every data point before making a single step (hence the name "batch"). An alternative approach to alleviate such computational costs is the Stochastic Gradient Descent method: essentially, the idea is to update the parameters $\theta$ sequentially (one data point at the time), with every observation $x_i, y_i$. Following the same notation we used for Gradient Descent, the following expression defines how to update parameters, while processing one data point at the time:$$ \theta^{(t+1)} = \theta^{(t)} - \eta \nabla_{\theta} J_i(\theta^{(t)}) $$The stochastic gradient approach allows us to start making progress on the minimization problem one step at the time. It is computationally cheaper, but it results in a larger variance of the loss function in comparison with batch gradient descent.Generally, the stochastic gradient descent method will get close to the optimal $\theta$ much faster than the batch method, but will never fully converge to the local (or global) minimum. Thus the stochastic gradient descent method is useful when we are satisfied with an **approximation** for the solution to our optimization problem. A full recipe for stochastic gradient descent follows:```python for i in range(n_epochs):  np.random.shuffle(data)  for example in data:    params_grad = evaluate_gradient(loss_function, example, params)    params = params - learning_rate * params_grad```The reshuffling of the data is done to avoid a bias in the optimization algorithm by providing the data examples in a particular order. Question 5. Implement your own version of Stochastic Gradient Descent, as a serial algorithm.Follow the guidelines below:     Define a function to perform gradient descent. The function should accept as inputs: the training data $x$ and $y$, the initial guess for the parameters ($\theta_1$ and $\theta_2$), the learning rate. Additional arguments include the definition of the maximum number of iterations before the algorithm stops, and a second stop condition on the marginal improvment on the loss.     Keep track of the values of the loss, for each iteration.     Keep track of the gradient values, for each iteration.Once the ```gradient_descent``` function is defined, you can generate input data according to the cell above, that use scikitlearn.    The output of your cell should contain the following information:     The values of the paramters obtained through Gradient Descent optimization      The values of the paramters obtained with the above cell, using scikitlearn      A plot of the loss versus iterations      A plot of the path the gradient takes from its initial to its final position<jupyter_code>def reShape(x):
    return np.c_[np.ones((x.shape[0])), x]


def loss(x,y,theta):
    lossRes=0
    for i in range (len(x)):
        lossRes+=(theta[0] + theta[1]*x[i]-y[i])**2
    return(lossRes/(2*len(x)))
    
def gradient(i,x,y,theta):
    res=x[i]*(theta[0]+theta[1]*x[i]-y[i])
    return(res)

def norm(t1,t2):
    return((t1[0]-t2[0])**2+(t1[1]-t2[1])**2)



def gradient_descent(x,Y,theta1, theta2,learning_rate, n_iteration, marginal_loss):
    X = reshape(x)
    theta=np.array([theta1, theta2])
    
    gradientRes = gradient(X,Y,theta)
    theta=theta.transpose()
    thetaNew=theta - learning_rate*gradientRes
    #print("theta",theta)
    
    lossRes=loss(x,Y,theta)
    
    lossHistory=[lossRes]
    gradientHistory=[]
    
    T1=[theta1]
    T2=[theta2]
    
    iteration=0
    

    
    while (iteration< n_iteration and norm(theta,thetaNew)>marginal_loss) : 
        
        for i in range(len(x)):     
            gradientRes = gradient(i,X,Y,theta)
            gradientHistory.append(gradientRes)    
            theta, thetaNew= thetaNew, thetaNew - learning_rate*gradientRes
            #print("theta","thetaNew",theta,thetaNew)      
            T1.append(theta[0])
            T2.append(theta[1])
            
            #print("loss",loss(x,Y,theta))
            lossHistory.append(loss(x,Y,thetaNew))
        
        iteration+=1
        
        
    return(lossHistory,gradientHistory,theta,T1,T2)


n_iteration=100
n_epochs=10

(lossHistory,gradientHistory,theta,T1,T2) = gradient_descent(x,y,1, 1,0.05, n_iteration, 0.0001)
  

print(theta)
print(type(lossHistory),type(gradientHistory))
print(len(lossHistory),len(gradientHistory))

abscisses=[i for i in range (len(lossHistory))]
print(type(abscisses))

plt.subplot(2, 1, 1)
plt.plot(abscisses, lossHistory)


print("abs: {} gradh:{}".format(abscisses, gradientHistory))
gradientX=[gradientHistory[i][0] for i in range(len(gradientHistory))]
gradientY=[gradientHistory[i][1] for i in range(len(gradientHistory))]
plt.subplot(2, 1, 2)

plt.plot(gradientX, gradientY)

plt.show()<jupyter_output>[305.55607711382902, 305.32153269380262, 305.21849147583771, 305.12581303507818, 305.04246737361461, 304.96751505424129, 304.9001106469907, 304.8394938715677, 304.78498104619297, 304.735957369652, 304.69186997993899, 304.65222171253163, 304.61656548779712, 304.58449926413118, 304.55566149976158, 304.52972707199098, 304.50640360766022, 304.48542818349819, 304.4665643589837, 304.44959950822516, 304.43434242072487, 304.42062114392354, 304.40828104313567, 304.39718305699017, 304.38720212864644, 304.37822579508929, 304.37015291856665, 304.36289254581936, 304.35636288227499, 304.35049036956696, 304.34520885599994, 304.34045885058555, 304.33618685222069, 304.33234474645025, 304.32888926299529, 304.32578148790856, 304.32298642490986, 304.3204726008575, 304.31821171099358, 304.31617829991046, 304.31434947463464] [-0.12511748  0.15056758]
[305.34518793847889, 305.14388429252415, 305.05549169758365, 304.97602269316809, 304.90458693490416, 304.84037266630793, 304.7826496774, 304.73076159089607, 30[...]
{'fraction_non_alphanumeric': 0.11154097372255439, 'fraction_numerical': 0.08260061297291865, 'mean_word_length': 4.453150684931507, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '223954', 'license': 'MIT', 'n_tokens_mistral': 7363, 'n_tokens_neox': 6180, 'n_words': 2133, 'path': 'Lab1-stochastic-gradient-descent/SGD-public.ipynb', 'repo_name': 'kanguyn/CLouds-Fall2017', 'size': '19903'}