diff --git a/doc/pub/week39/html/week39-bs.html b/doc/pub/week39/html/week39-bs.html index 165880adb..4d8240a34 100644 --- a/doc/pub/week39/html/week39-bs.html +++ b/doc/pub/week39/html/week39-bs.html @@ -236,6 +236,10 @@ doconce format html week39.do.txt --html_style=bootstrap --pygments_html_style=d 2, None, 'same-code-but-now-with-momentum-gradient-descent'), + ("But noen of these can compete with Newton's method", + 2, + None, + 'but-noen-of-these-can-compete-with-newton-s-method'), ('Including Stochastic Gradient Descent with Autograd', 2, None, @@ -244,11 +248,16 @@ doconce format html week39.do.txt --html_style=bootstrap --pygments_html_style=d 2, None, 'same-code-but-now-with-momentum-gradient-descent'), + ('Same problem but now with AdaGrad', + 2, + None, + 'same-problem-but-now-with-adagrad'), ('And Logistic Regression', 2, None, 'and-logistic-regression'), ('Introducing "JAX":"https://jax.readthedocs.io/en/latest/"', 2, None, - 'introducing-jax-https-jax-readthedocs-io-en-latest')]} + 'introducing-jax-https-jax-readthedocs-io-en-latest'), + ('Weekend challenge', 2, None, 'weekend-challenge')]} end of tocinfo --> @@ -327,7 +336,7 @@ MathJax.Hub.Config({
  • Program example for gradient descent with Ridge Regression
  • Using gradient descent methods, limitations
  • Improving gradient descent with momentum
  • -
  • Same code but now with momentum gradient descent
  • +
  • Same code but now with momentum gradient descent
  • Overview video on Stochastic Gradient Descent
  • Batches and mini-batches
  • Stochastic Gradient Descent (SGD)
  • @@ -360,11 +369,14 @@ MathJax.Hub.Config({
  • The syntax a.dot(b) when finding the dot product
  • Recommended to avoid
  • Using Autograd with OLS
  • -
  • Same code but now with momentum gradient descent
  • -
  • Including Stochastic Gradient Descent with Autograd
  • -
  • Same code but now with momentum gradient descent
  • -
  • And Logistic Regression
  • -
  • Introducing "JAX":"https://jax.readthedocs.io/en/latest/"
  • +
  • Same code but now with momentum gradient descent
  • +
  • But noen of these can compete with Newton's method
  • +
  • Including Stochastic Gradient Descent with Autograd
  • +
  • Same code but now with momentum gradient descent
  • +
  • Same problem but now with AdaGrad
  • +
  • And Logistic Regression
  • +
  • Introducing "JAX":"https://jax.readthedocs.io/en/latest/"
  • +
  • Weekend challenge
  • @@ -419,7 +431,7 @@ MathJax.Hub.Config({
  • 9
  • 10
  • ...
  • -
  • 83
  • +
  • 86
  • »
  • diff --git a/doc/pub/week39/html/week39-reveal.html b/doc/pub/week39/html/week39-reveal.html index 0a0f5455a..31fd0418e 100644 --- a/doc/pub/week39/html/week39-reveal.html +++ b/doc/pub/week39/html/week39-reveal.html @@ -3229,6 +3229,70 @@ delta_momentum = 0.3 +
    +

    But noen of these can compete with Newton's method

    + + + +
    +
    +
    +
    +
    +
    # Using Newton's method
    +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
    +beta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(beta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +# Note that here the Hessian does not depend on the parameters beta
    +invH = np.linalg.pinv(H)
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +beta = np.random.randn(2,1)
    +Niterations = 5
    +
    +# define the gradient
    +training_gradient = grad(CostOLS)
    +
    +for iter in range(Niterations):
    +    gradients = training_gradient(beta)
    +    beta -= invH @ gradients
    +    print(iter,gradients[0],gradients[1])
    +print("beta from own Newton code")
    +print(beta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +

    Including Stochastic Gradient Descent with Autograd

    In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using autograd.

    @@ -3421,6 +3485,92 @@ delta_momentum = 0.3
    +
    +

    Same problem but now with AdaGrad

    + + +
    +
    +
    +
    +
    +
    # 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
    +
    +# Note change from previous example
    +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 = 100
    +
    +# 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)
    +print(np.size(gradients))
    +
    +n_epochs = 50
    +M = 5   #size of each minibatch
    +m = int(n/M) #number of minibatches
    +theta = np.random.randn(2,1)
    +
    +
    +# Including AdaGrad
    +delta  = 0.000001
    +for epoch in range(n_epochs):
    +    for i in range(m):
    +        random_index = M*np.random.randint(m)
    +        xi = X[random_index:random_index+M]
    +        yi = y[random_index:random_index+M]
    +        gradients = (1.0/M)*training_gradient(yi, xi, theta)
    +        # calculate squared gradient by Hadamard multiplication
    +        r -= gradients*gradients
    +        # compute update
    +        update = (1.0/delta+np.sqrt(r))*gradients
    +        theta = eta*update
    +print("theta from own AdaGrad")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +

    And Logistic Regression

    @@ -3524,6 +3674,16 @@ derivative_fn = grad(sum_logistic)
    +
    +

    Weekend challenge

    + + +
    + diff --git a/doc/pub/week39/html/week39-solarized.html b/doc/pub/week39/html/week39-solarized.html index 0cdc418fa..76bc2d149 100644 --- a/doc/pub/week39/html/week39-solarized.html +++ b/doc/pub/week39/html/week39-solarized.html @@ -263,6 +263,10 @@ div.toc p,a { 2, None, 'same-code-but-now-with-momentum-gradient-descent'), + ("But noen of these can compete with Newton's method", + 2, + None, + 'but-noen-of-these-can-compete-with-newton-s-method'), ('Including Stochastic Gradient Descent with Autograd', 2, None, @@ -271,11 +275,16 @@ div.toc p,a { 2, None, 'same-code-but-now-with-momentum-gradient-descent'), + ('Same problem but now with AdaGrad', + 2, + None, + 'same-problem-but-now-with-adagrad'), ('And Logistic Regression', 2, None, 'and-logistic-regression'), ('Introducing "JAX":"https://jax.readthedocs.io/en/latest/"', 2, None, - 'introducing-jax-https-jax-readthedocs-io-en-latest')]} + 'introducing-jax-https-jax-readthedocs-io-en-latest'), + ('Weekend challenge', 2, None, 'weekend-challenge')]} end of tocinfo --> @@ -3147,6 +3156,70 @@ delta_momentum = 0.3 +









    +

    But noen of these can compete with Newton's method

    + + + +
    +
    +
    +
    +
    +
    # Using Newton's method
    +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
    +beta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(beta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +# Note that here the Hessian does not depend on the parameters beta
    +invH = np.linalg.pinv(H)
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +beta = np.random.randn(2,1)
    +Niterations = 5
    +
    +# define the gradient
    +training_gradient = grad(CostOLS)
    +
    +for iter in range(Niterations):
    +    gradients = training_gradient(beta)
    +    beta -= invH @ gradients
    +    print(iter,gradients[0],gradients[1])
    +print("beta from own Newton code")
    +print(beta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + +









    Including Stochastic Gradient Descent with Autograd

    In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using autograd.

    @@ -3339,6 +3412,92 @@ delta_momentum = 0.3 +









    +

    Same problem but now with AdaGrad

    + + +
    +
    +
    +
    +
    +
    # 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
    +
    +# Note change from previous example
    +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 = 100
    +
    +# 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)
    +print(np.size(gradients))
    +
    +n_epochs = 50
    +M = 5   #size of each minibatch
    +m = int(n/M) #number of minibatches
    +theta = np.random.randn(2,1)
    +
    +
    +# Including AdaGrad
    +delta  = 0.000001
    +for epoch in range(n_epochs):
    +    for i in range(m):
    +        random_index = M*np.random.randint(m)
    +        xi = X[random_index:random_index+M]
    +        yi = y[random_index:random_index+M]
    +        gradients = (1.0/M)*training_gradient(yi, xi, theta)
    +        # calculate squared gradient by Hadamard multiplication
    +        r -= gradients*gradients
    +        # compute update
    +        update = (1.0/delta+np.sqrt(r))*gradients
    +        theta = eta*update
    +print("theta from own AdaGrad")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + +









    And Logistic Regression

    @@ -3441,6 +3600,15 @@ derivative_fn = grad(sum_logistic) + +









    +

    Weekend challenge

    + +
    © 1999-2022, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license diff --git a/doc/pub/week39/html/week39.html b/doc/pub/week39/html/week39.html index e2a59f217..953335216 100644 --- a/doc/pub/week39/html/week39.html +++ b/doc/pub/week39/html/week39.html @@ -340,6 +340,10 @@ div.toc p,a { 2, None, 'same-code-but-now-with-momentum-gradient-descent'), + ("But noen of these can compete with Newton's method", + 2, + None, + 'but-noen-of-these-can-compete-with-newton-s-method'), ('Including Stochastic Gradient Descent with Autograd', 2, None, @@ -348,11 +352,16 @@ div.toc p,a { 2, None, 'same-code-but-now-with-momentum-gradient-descent'), + ('Same problem but now with AdaGrad', + 2, + None, + 'same-problem-but-now-with-adagrad'), ('And Logistic Regression', 2, None, 'and-logistic-regression'), ('Introducing "JAX":"https://jax.readthedocs.io/en/latest/"', 2, None, - 'introducing-jax-https-jax-readthedocs-io-en-latest')]} + 'introducing-jax-https-jax-readthedocs-io-en-latest'), + ('Weekend challenge', 2, None, 'weekend-challenge')]} end of tocinfo --> @@ -3224,6 +3233,70 @@ delta_momentum = But noen of these can compete with Newton's method + + + +
    +
    +
    +
    +
    +
    # Using Newton's method
    +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
    +beta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(beta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +# Note that here the Hessian does not depend on the parameters beta
    +invH = np.linalg.pinv(H)
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +beta = np.random.randn(2,1)
    +Niterations = 5
    +
    +# define the gradient
    +training_gradient = grad(CostOLS)
    +
    +for iter in range(Niterations):
    +    gradients = training_gradient(beta)
    +    beta -= invH @ gradients
    +    print(iter,gradients[0],gradients[1])
    +print("beta from own Newton code")
    +print(beta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + +









    Including Stochastic Gradient Descent with Autograd

    In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using autograd.

    @@ -3416,6 +3489,92 @@ delta_momentum = Same problem but now with AdaGrad + + +
    +
    +
    +
    +
    +
    # 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
    +
    +# Note change from previous example
    +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 = 100
    +
    +# 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)
    +print(np.size(gradients))
    +
    +n_epochs = 50
    +M = 5   #size of each minibatch
    +m = int(n/M) #number of minibatches
    +theta = np.random.randn(2,1)
    +
    +
    +# Including AdaGrad
    +delta  = 0.000001
    +for epoch in range(n_epochs):
    +    for i in range(m):
    +        random_index = M*np.random.randint(m)
    +        xi = X[random_index:random_index+M]
    +        yi = y[random_index:random_index+M]
    +        gradients = (1.0/M)*training_gradient(yi, xi, theta)
    +        # calculate squared gradient by Hadamard multiplication
    +        r -= gradients*gradients
    +        # compute update
    +        update = (1.0/delta+np.sqrt(r))*gradients
    +        theta = eta*update
    +print("theta from own AdaGrad")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + +









    And Logistic Regression

    @@ -3518,6 +3677,15 @@ derivative_fn = grad(sum_logistic) + +









    +

    Weekend challenge

    + +
      +
    • Try to run the above codes and implement the stochastic gradient descent with RMSprop and ADAM.
    • +
    • Add a more complicated function and study the rate of convergence for the derivatives as function of the different methods
    • +
    • Extend from linear regression to logistic regression.
    • +
    © 1999-2022, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license diff --git a/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz b/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz index 3b8729668..9b9d1c2cb 100644 Binary files a/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz and b/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz differ diff --git a/doc/pub/week39/ipynb/week39.ipynb b/doc/pub/week39/ipynb/week39.ipynb index 181035900..1e8f64dce 100644 --- a/doc/pub/week39/ipynb/week39.ipynb +++ b/doc/pub/week39/ipynb/week39.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "a3f2a49b", + "id": "3090dab9", "metadata": { "editable": true }, @@ -14,7 +14,7 @@ }, { "cell_type": "markdown", - "id": "ad1e01ff", + "id": "ab0dea94", "metadata": { "editable": true }, @@ -29,7 +29,7 @@ }, { "cell_type": "markdown", - "id": "1ce0a810", + "id": "9db200b9", "metadata": { "editable": true }, @@ -55,7 +55,7 @@ }, { "cell_type": "markdown", - "id": "8ea463cd", + "id": "daa1efbf", "metadata": { "editable": true }, @@ -76,7 +76,7 @@ }, { "cell_type": "markdown", - "id": "7bcc3cee", + "id": "d67302c7", "metadata": { "editable": true }, @@ -93,7 +93,7 @@ }, { "cell_type": "markdown", - "id": "8258a91e", + "id": "6cffc694", "metadata": { "editable": true }, @@ -108,7 +108,7 @@ }, { "cell_type": "markdown", - "id": "fc4d239c", + "id": "7c0e4b41", "metadata": { "editable": true }, @@ -118,7 +118,7 @@ }, { "cell_type": "markdown", - "id": "5e411ed9", + "id": "a4d042c9", "metadata": { "editable": true }, @@ -134,7 +134,7 @@ }, { "cell_type": "markdown", - "id": "8559bad7", + "id": "a2e6e6ca", "metadata": { "editable": true }, @@ -146,7 +146,7 @@ }, { "cell_type": "markdown", - "id": "8cb89fe1", + "id": "44ba6d7f", "metadata": { "editable": true }, @@ -157,7 +157,7 @@ }, { "cell_type": "markdown", - "id": "4e7f874e", + "id": "f83dfe9e", "metadata": { "editable": true }, @@ -169,7 +169,7 @@ }, { "cell_type": "markdown", - "id": "1f6c709e", + "id": "36ae9e1f", "metadata": { "editable": true }, @@ -179,7 +179,7 @@ }, { "cell_type": "markdown", - "id": "618e2feb", + "id": "f68dff87", "metadata": { "editable": true }, @@ -193,7 +193,7 @@ }, { "cell_type": "markdown", - "id": "19b17dae", + "id": "1c04ada4", "metadata": { "editable": true }, @@ -205,7 +205,7 @@ }, { "cell_type": "markdown", - "id": "9ad58b9c", + "id": "57d59ede", "metadata": { "editable": true }, @@ -215,7 +215,7 @@ }, { "cell_type": "markdown", - "id": "6763d42d", + "id": "6731c775", "metadata": { "editable": true }, @@ -227,7 +227,7 @@ }, { "cell_type": "markdown", - "id": "ce95b330", + "id": "344395d0", "metadata": { "editable": true }, @@ -239,7 +239,7 @@ }, { "cell_type": "markdown", - "id": "4c7a5561", + "id": "030a7dd9", "metadata": { "editable": true }, @@ -259,7 +259,7 @@ }, { "cell_type": "markdown", - "id": "b2f5f204", + "id": "f714706e", "metadata": { "editable": true }, @@ -275,7 +275,7 @@ }, { "cell_type": "markdown", - "id": "25f1ec55", + "id": "bc0a90df", "metadata": { "editable": true }, @@ -291,7 +291,7 @@ }, { "cell_type": "markdown", - "id": "12800af1", + "id": "2ced70dc", "metadata": { "editable": true }, @@ -302,7 +302,7 @@ }, { "cell_type": "markdown", - "id": "24eac4c2", + "id": "f69cd4d8", "metadata": { "editable": true }, @@ -314,7 +314,7 @@ }, { "cell_type": "markdown", - "id": "fefe799b", + "id": "9478068d", "metadata": { "editable": true }, @@ -324,7 +324,7 @@ }, { "cell_type": "markdown", - "id": "ac788ad2", + "id": "66e95761", "metadata": { "editable": true }, @@ -336,7 +336,7 @@ }, { "cell_type": "markdown", - "id": "6837c77b", + "id": "3b1f9fb4", "metadata": { "editable": true }, @@ -346,7 +346,7 @@ }, { "cell_type": "markdown", - "id": "ea620a33", + "id": "26695601", "metadata": { "editable": true }, @@ -358,7 +358,7 @@ }, { "cell_type": "markdown", - "id": "12439301", + "id": "1951019b", "metadata": { "editable": true }, @@ -380,7 +380,7 @@ }, { "cell_type": "markdown", - "id": "2991e70c", + "id": "b38f86e5", "metadata": { "editable": true }, @@ -393,7 +393,7 @@ }, { "cell_type": "markdown", - "id": "0bb42d16", + "id": "81a4e8bd", "metadata": { "editable": true }, @@ -406,7 +406,7 @@ }, { "cell_type": "markdown", - "id": "cc166498", + "id": "12f207a8", "metadata": { "editable": true }, @@ -416,7 +416,7 @@ }, { "cell_type": "markdown", - "id": "941aae4f", + "id": "440c1a73", "metadata": { "editable": true }, @@ -434,7 +434,7 @@ }, { "cell_type": "markdown", - "id": "957408a1", + "id": "9e301758", "metadata": { "editable": true }, @@ -444,7 +444,7 @@ }, { "cell_type": "markdown", - "id": "3d6b3d4a", + "id": "714caf32", "metadata": { "editable": true }, @@ -459,7 +459,7 @@ }, { "cell_type": "markdown", - "id": "17d6443c", + "id": "79945812", "metadata": { "editable": true }, @@ -469,7 +469,7 @@ }, { "cell_type": "markdown", - "id": "583a3f9c", + "id": "821bfe43", "metadata": { "editable": true }, @@ -483,7 +483,7 @@ }, { "cell_type": "markdown", - "id": "55065669", + "id": "7a329d62", "metadata": { "editable": true }, @@ -493,7 +493,7 @@ }, { "cell_type": "markdown", - "id": "8ca8e144", + "id": "a759ca7a", "metadata": { "editable": true }, @@ -507,7 +507,7 @@ }, { "cell_type": "markdown", - "id": "9173c236", + "id": "37bd514e", "metadata": { "editable": true }, @@ -522,7 +522,7 @@ }, { "cell_type": "markdown", - "id": "34850c9f", + "id": "aa97f992", "metadata": { "editable": true }, @@ -539,7 +539,7 @@ }, { "cell_type": "markdown", - "id": "524b6a67", + "id": "67fa21b3", "metadata": { "editable": true }, @@ -551,7 +551,7 @@ }, { "cell_type": "markdown", - "id": "b1d69a19", + "id": "f8f33cfa", "metadata": { "editable": true }, @@ -565,7 +565,7 @@ }, { "cell_type": "markdown", - "id": "5c48c377", + "id": "12fcbcb4", "metadata": { "editable": true }, @@ -580,7 +580,7 @@ }, { "cell_type": "markdown", - "id": "d847b012", + "id": "247f240a", "metadata": { "editable": true }, @@ -592,7 +592,7 @@ }, { "cell_type": "markdown", - "id": "430c26cb", + "id": "a7f4feb0", "metadata": { "editable": true }, @@ -603,7 +603,7 @@ }, { "cell_type": "markdown", - "id": "988b26d1", + "id": "581dc6c7", "metadata": { "editable": true }, @@ -631,7 +631,7 @@ }, { "cell_type": "markdown", - "id": "5527bcc6", + "id": "115cc2f5", "metadata": { "editable": true }, @@ -653,7 +653,7 @@ }, { "cell_type": "markdown", - "id": "85a3c5fc", + "id": "981820fd", "metadata": { "editable": true }, @@ -675,7 +675,7 @@ }, { "cell_type": "markdown", - "id": "103ca9e0", + "id": "71769268", "metadata": { "editable": true }, @@ -687,7 +687,7 @@ }, { "cell_type": "markdown", - "id": "5f043b19", + "id": "7236d516", "metadata": { "editable": true }, @@ -724,7 +724,7 @@ }, { "cell_type": "markdown", - "id": "7c5ac593", + "id": "1be5600b", "metadata": { "editable": true }, @@ -752,7 +752,7 @@ }, { "cell_type": "markdown", - "id": "6bd0a3f3", + "id": "bc9479dd", "metadata": { "editable": true }, @@ -782,7 +782,7 @@ }, { "cell_type": "markdown", - "id": "3e40ce1f", + "id": "b287dc6b", "metadata": { "editable": true }, @@ -802,7 +802,7 @@ }, { "cell_type": "markdown", - "id": "ade236c4", + "id": "f3b051c9", "metadata": { "editable": true }, @@ -814,7 +814,7 @@ }, { "cell_type": "markdown", - "id": "e1c7967b", + "id": "1d60b0bb", "metadata": { "editable": true }, @@ -824,7 +824,7 @@ }, { "cell_type": "markdown", - "id": "ad63044a", + "id": "cee9ebc6", "metadata": { "editable": true }, @@ -836,7 +836,7 @@ }, { "cell_type": "markdown", - "id": "2151a6e5", + "id": "c069bba3", "metadata": { "editable": true }, @@ -848,7 +848,7 @@ }, { "cell_type": "markdown", - "id": "08d1dcfc", + "id": "7c83fe7f", "metadata": { "editable": true }, @@ -860,7 +860,7 @@ }, { "cell_type": "markdown", - "id": "714a13d7", + "id": "cda7891d", "metadata": { "editable": true }, @@ -872,7 +872,7 @@ }, { "cell_type": "markdown", - "id": "5447991e", + "id": "7edb97a2", "metadata": { "editable": true }, @@ -883,7 +883,7 @@ }, { "cell_type": "markdown", - "id": "9d8d5636", + "id": "aab86b4c", "metadata": { "editable": true }, @@ -896,7 +896,7 @@ }, { "cell_type": "markdown", - "id": "c4c62bf8", + "id": "6d5fa655", "metadata": { "editable": true }, @@ -908,7 +908,7 @@ }, { "cell_type": "markdown", - "id": "e75ef87c", + "id": "551ffdae", "metadata": { "editable": true }, @@ -918,7 +918,7 @@ }, { "cell_type": "markdown", - "id": "967df7e5", + "id": "697ac0c7", "metadata": { "editable": true }, @@ -930,7 +930,7 @@ }, { "cell_type": "markdown", - "id": "9770216b", + "id": "2501fd5b", "metadata": { "editable": true }, @@ -940,7 +940,7 @@ }, { "cell_type": "markdown", - "id": "600bb9bf", + "id": "3ef8cf2e", "metadata": { "editable": true }, @@ -951,7 +951,7 @@ }, { "cell_type": "markdown", - "id": "1172ae79", + "id": "2cb81217", "metadata": { "editable": true }, @@ -963,7 +963,7 @@ }, { "cell_type": "markdown", - "id": "3760450e", + "id": "3f57e6c0", "metadata": { "editable": true }, @@ -975,7 +975,7 @@ }, { "cell_type": "markdown", - "id": "b0441a59", + "id": "e4f4feb2", "metadata": { "editable": true }, @@ -987,7 +987,7 @@ }, { "cell_type": "markdown", - "id": "88c2cdab", + "id": "e311ea5a", "metadata": { "editable": true }, @@ -998,7 +998,7 @@ }, { "cell_type": "markdown", - "id": "df8ddc4b", + "id": "ccdf1d1f", "metadata": { "editable": true }, @@ -1009,7 +1009,7 @@ }, { "cell_type": "markdown", - "id": "0f297c3d", + "id": "032554f8", "metadata": { "editable": true }, @@ -1021,7 +1021,7 @@ }, { "cell_type": "markdown", - "id": "0d0d0aa2", + "id": "e5278ef0", "metadata": { "editable": true }, @@ -1031,7 +1031,7 @@ }, { "cell_type": "markdown", - "id": "c4d38149", + "id": "0abff66a", "metadata": { "editable": true }, @@ -1043,7 +1043,7 @@ }, { "cell_type": "markdown", - "id": "d4372560", + "id": "aaa20c9f", "metadata": { "editable": true }, @@ -1053,7 +1053,7 @@ }, { "cell_type": "markdown", - "id": "e718ad83", + "id": "88de1b2c", "metadata": { "editable": true }, @@ -1065,7 +1065,7 @@ }, { "cell_type": "markdown", - "id": "457ea769", + "id": "a6b7e7eb", "metadata": { "editable": true }, @@ -1075,7 +1075,7 @@ }, { "cell_type": "markdown", - "id": "e5677842", + "id": "f48a1322", "metadata": { "editable": true }, @@ -1087,7 +1087,7 @@ }, { "cell_type": "markdown", - "id": "f41dcf2f", + "id": "fde12a28", "metadata": { "editable": true }, @@ -1097,7 +1097,7 @@ }, { "cell_type": "markdown", - "id": "c3a84d60", + "id": "47d4f0f8", "metadata": { "editable": true }, @@ -1109,7 +1109,7 @@ }, { "cell_type": "markdown", - "id": "54eec9d8", + "id": "1571831b", "metadata": { "editable": true }, @@ -1120,7 +1120,7 @@ { "cell_type": "code", "execution_count": 1, - "id": "e65d59f3", + "id": "dcc70cc9", "metadata": { "collapsed": false, "editable": true @@ -1153,7 +1153,7 @@ }, { "cell_type": "markdown", - "id": "a1b1ec51", + "id": "91654b7a", "metadata": { "editable": true }, @@ -1164,7 +1164,7 @@ { "cell_type": "code", "execution_count": 2, - "id": "d509c489", + "id": "43f2a763", "metadata": { "collapsed": false, "editable": true @@ -1178,7 +1178,7 @@ }, { "cell_type": "markdown", - "id": "0a234ac1", + "id": "f6ad5e15", "metadata": { "editable": true }, @@ -1189,7 +1189,7 @@ { "cell_type": "code", "execution_count": 3, - "id": "6dee8ee2", + "id": "be373c58", "metadata": { "collapsed": false, "editable": true @@ -1202,7 +1202,7 @@ }, { "cell_type": "markdown", - "id": "6739f9d4", + "id": "e48b0f67", "metadata": { "editable": true }, @@ -1213,7 +1213,7 @@ { "cell_type": "code", "execution_count": 4, - "id": "a635a93f", + "id": "bc6ae2a9", "metadata": { "collapsed": false, "editable": true @@ -1231,7 +1231,7 @@ }, { "cell_type": "markdown", - "id": "1738466f", + "id": "5618d37d", "metadata": { "editable": true }, @@ -1242,7 +1242,7 @@ { "cell_type": "code", "execution_count": 5, - "id": "c58d0115", + "id": "0d22b148", "metadata": { "collapsed": false, "editable": true @@ -1257,7 +1257,7 @@ }, { "cell_type": "markdown", - "id": "7ea54653", + "id": "98c2a710", "metadata": { "editable": true }, @@ -1267,7 +1267,7 @@ }, { "cell_type": "markdown", - "id": "dd7e00ed", + "id": "2b589600", "metadata": { "editable": true }, @@ -1281,7 +1281,7 @@ }, { "cell_type": "markdown", - "id": "7f01fc09", + "id": "8cb6ed9e", "metadata": { "editable": true }, @@ -1293,7 +1293,7 @@ }, { "cell_type": "markdown", - "id": "6205a0c0", + "id": "d05dc72f", "metadata": { "editable": true }, @@ -1304,7 +1304,7 @@ }, { "cell_type": "markdown", - "id": "0a474db3", + "id": "612657ad", "metadata": { "editable": true }, @@ -1316,7 +1316,7 @@ }, { "cell_type": "markdown", - "id": "281c3e87", + "id": "c14e7863", "metadata": { "editable": true }, @@ -1327,7 +1327,7 @@ }, { "cell_type": "markdown", - "id": "72adadf6", + "id": "54f184de", "metadata": { "editable": true }, @@ -1338,7 +1338,7 @@ }, { "cell_type": "markdown", - "id": "b1a60997", + "id": "acd6d248", "metadata": { "editable": true }, @@ -1350,7 +1350,7 @@ }, { "cell_type": "markdown", - "id": "4d655399", + "id": "24083680", "metadata": { "editable": true }, @@ -1360,7 +1360,7 @@ }, { "cell_type": "markdown", - "id": "116ebbac", + "id": "b41beb62", "metadata": { "editable": true }, @@ -1372,7 +1372,7 @@ }, { "cell_type": "markdown", - "id": "55185d8f", + "id": "2626b306", "metadata": { "editable": true }, @@ -1384,7 +1384,7 @@ }, { "cell_type": "markdown", - "id": "afd187d3", + "id": "9c485499", "metadata": { "editable": true }, @@ -1396,7 +1396,7 @@ }, { "cell_type": "markdown", - "id": "06b67a8c", + "id": "d35c34d6", "metadata": { "editable": true }, @@ -1408,7 +1408,7 @@ }, { "cell_type": "markdown", - "id": "a6c814bc", + "id": "4ed23f34", "metadata": { "editable": true }, @@ -1419,7 +1419,7 @@ }, { "cell_type": "markdown", - "id": "8e4c87d2", + "id": "f9a65949", "metadata": { "editable": true }, @@ -1431,7 +1431,7 @@ }, { "cell_type": "markdown", - "id": "71ac6a87", + "id": "41497653", "metadata": { "editable": true }, @@ -1441,7 +1441,7 @@ }, { "cell_type": "markdown", - "id": "2a5a1639", + "id": "218bf4a2", "metadata": { "editable": true }, @@ -1453,7 +1453,7 @@ }, { "cell_type": "markdown", - "id": "a21020d2", + "id": "f903700f", "metadata": { "editable": true }, @@ -1463,7 +1463,7 @@ }, { "cell_type": "markdown", - "id": "fc5b548c", + "id": "05bf8370", "metadata": { "editable": true }, @@ -1475,7 +1475,7 @@ }, { "cell_type": "markdown", - "id": "52b40bde", + "id": "afea8ad1", "metadata": { "editable": true }, @@ -1495,7 +1495,7 @@ }, { "cell_type": "markdown", - "id": "361660ca", + "id": "0605b9e9", "metadata": { "editable": true }, @@ -1507,7 +1507,7 @@ }, { "cell_type": "markdown", - "id": "519d134c", + "id": "ca8bc8e7", "metadata": { "editable": true }, @@ -1517,7 +1517,7 @@ }, { "cell_type": "markdown", - "id": "0a7faed4", + "id": "2a752ed8", "metadata": { "editable": true }, @@ -1529,7 +1529,7 @@ }, { "cell_type": "markdown", - "id": "e7a1c68e", + "id": "64932658", "metadata": { "editable": true }, @@ -1539,7 +1539,7 @@ }, { "cell_type": "markdown", - "id": "6226cc6a", + "id": "1e8db58a", "metadata": { "editable": true }, @@ -1550,7 +1550,7 @@ }, { "cell_type": "markdown", - "id": "2e7286a9", + "id": "417388f7", "metadata": { "editable": true }, @@ -1562,7 +1562,7 @@ }, { "cell_type": "markdown", - "id": "356d6972", + "id": "a2884ba6", "metadata": { "editable": true }, @@ -1574,7 +1574,7 @@ }, { "cell_type": "markdown", - "id": "b8d473e6", + "id": "bafbb7bc", "metadata": { "editable": true }, @@ -1586,7 +1586,7 @@ }, { "cell_type": "markdown", - "id": "2d22f619", + "id": "335831d3", "metadata": { "editable": true }, @@ -1599,7 +1599,7 @@ }, { "cell_type": "markdown", - "id": "16102292", + "id": "ca8f091c", "metadata": { "editable": true }, @@ -1610,7 +1610,7 @@ }, { "cell_type": "markdown", - "id": "8cac320f", + "id": "d353c7de", "metadata": { "editable": true }, @@ -1622,7 +1622,7 @@ }, { "cell_type": "markdown", - "id": "7d17a555", + "id": "97ad2759", "metadata": { "editable": true }, @@ -1638,7 +1638,7 @@ }, { "cell_type": "markdown", - "id": "dfc11703", + "id": "8e37f2a6", "metadata": { "editable": true }, @@ -1650,7 +1650,7 @@ }, { "cell_type": "markdown", - "id": "511d4f6c", + "id": "988d5180", "metadata": { "editable": true }, @@ -1661,7 +1661,7 @@ }, { "cell_type": "markdown", - "id": "167ad0dd", + "id": "40dddf68", "metadata": { "editable": true }, @@ -1673,7 +1673,7 @@ }, { "cell_type": "markdown", - "id": "5e3e8f17", + "id": "6e30b8e8", "metadata": { "editable": true }, @@ -1683,7 +1683,7 @@ }, { "cell_type": "markdown", - "id": "d17a807f", + "id": "6f9375f6", "metadata": { "editable": true }, @@ -1695,7 +1695,7 @@ }, { "cell_type": "markdown", - "id": "c6545761", + "id": "429f93b7", "metadata": { "editable": true }, @@ -1705,7 +1705,7 @@ }, { "cell_type": "markdown", - "id": "48fa9c57", + "id": "50a29c3f", "metadata": { "editable": true }, @@ -1717,7 +1717,7 @@ }, { "cell_type": "markdown", - "id": "14ca3319", + "id": "54b74e4b", "metadata": { "editable": true }, @@ -1727,7 +1727,7 @@ }, { "cell_type": "markdown", - "id": "96480b6a", + "id": "22278fd2", "metadata": { "editable": true }, @@ -1739,7 +1739,7 @@ }, { "cell_type": "markdown", - "id": "f5f3cbcf", + "id": "4ef8149e", "metadata": { "editable": true }, @@ -1763,7 +1763,7 @@ { "cell_type": "code", "execution_count": 6, - "id": "58aecda1", + "id": "9ab0e233", "metadata": { "collapsed": false, "editable": true @@ -1776,7 +1776,7 @@ }, { "cell_type": "markdown", - "id": "307bbb57", + "id": "c4160ae4", "metadata": { "editable": true }, @@ -1787,7 +1787,7 @@ }, { "cell_type": "markdown", - "id": "7a110c79", + "id": "89a9d484", "metadata": { "editable": true }, @@ -1799,7 +1799,7 @@ }, { "cell_type": "markdown", - "id": "28b08d11", + "id": "bca24c47", "metadata": { "editable": true }, @@ -1809,7 +1809,7 @@ }, { "cell_type": "markdown", - "id": "c28cd1d6", + "id": "9218a2d9", "metadata": { "editable": true }, @@ -1821,7 +1821,7 @@ }, { "cell_type": "markdown", - "id": "4e123dce", + "id": "527cbe33", "metadata": { "editable": true }, @@ -1835,7 +1835,7 @@ }, { "cell_type": "markdown", - "id": "d4dae5a6", + "id": "a98916d2", "metadata": { "editable": true }, @@ -1851,7 +1851,7 @@ }, { "cell_type": "markdown", - "id": "f02bdc7a", + "id": "511077bb", "metadata": { "editable": true }, @@ -1861,7 +1861,7 @@ }, { "cell_type": "markdown", - "id": "8a27caf3", + "id": "93fcdfde", "metadata": { "editable": true }, @@ -1873,7 +1873,7 @@ }, { "cell_type": "markdown", - "id": "750e00b1", + "id": "0e01b869", "metadata": { "editable": true }, @@ -1883,7 +1883,7 @@ }, { "cell_type": "markdown", - "id": "24b9613d", + "id": "881fedf9", "metadata": { "editable": true }, @@ -1895,7 +1895,7 @@ }, { "cell_type": "markdown", - "id": "68208c7e", + "id": "3fd08058", "metadata": { "editable": true }, @@ -1909,7 +1909,7 @@ }, { "cell_type": "markdown", - "id": "456d3441", + "id": "90d9a9f3", "metadata": { "editable": true }, @@ -1919,7 +1919,7 @@ }, { "cell_type": "markdown", - "id": "5e9d90e1", + "id": "f52a1b2f", "metadata": { "editable": true }, @@ -1930,7 +1930,7 @@ }, { "cell_type": "markdown", - "id": "99072c22", + "id": "9c636475", "metadata": { "editable": true }, @@ -1945,7 +1945,7 @@ }, { "cell_type": "markdown", - "id": "339105e4", + "id": "1efbe493", "metadata": { "editable": true }, @@ -1955,7 +1955,7 @@ }, { "cell_type": "markdown", - "id": "fa3aeb52", + "id": "6f68a417", "metadata": { "editable": true }, @@ -1967,7 +1967,7 @@ }, { "cell_type": "markdown", - "id": "868d976a", + "id": "ceed79ec", "metadata": { "editable": true }, @@ -1979,7 +1979,7 @@ }, { "cell_type": "markdown", - "id": "34207df7", + "id": "9a16d440", "metadata": { "editable": true }, @@ -1994,7 +1994,7 @@ }, { "cell_type": "markdown", - "id": "38a9d087", + "id": "1d6881a7", "metadata": { "editable": true }, @@ -2007,7 +2007,7 @@ { "cell_type": "code", "execution_count": 7, - "id": "dc7badb6", + "id": "70a26a35", "metadata": { "collapsed": false, "editable": true @@ -2064,7 +2064,7 @@ }, { "cell_type": "markdown", - "id": "59d83633", + "id": "f6e097fa", "metadata": { "editable": true }, @@ -2075,7 +2075,7 @@ { "cell_type": "code", "execution_count": 8, - "id": "ccaa3e15", + "id": "5f932e9e", "metadata": { "collapsed": false, "editable": true @@ -2102,7 +2102,7 @@ }, { "cell_type": "markdown", - "id": "b6683bd2", + "id": "23a2fcd9", "metadata": { "editable": true }, @@ -2114,7 +2114,7 @@ }, { "cell_type": "markdown", - "id": "2bb6c0fe", + "id": "37e42d79", "metadata": { "editable": true }, @@ -2126,7 +2126,7 @@ }, { "cell_type": "markdown", - "id": "e9d915fa", + "id": "6f372e4c", "metadata": { "editable": true }, @@ -2136,7 +2136,7 @@ }, { "cell_type": "markdown", - "id": "f05f103f", + "id": "a6216267", "metadata": { "editable": true }, @@ -2150,7 +2150,7 @@ }, { "cell_type": "markdown", - "id": "f1006d59", + "id": "97566e5a", "metadata": { "editable": true }, @@ -2160,7 +2160,7 @@ }, { "cell_type": "markdown", - "id": "64482546", + "id": "c4ddee5f", "metadata": { "editable": true }, @@ -2172,7 +2172,7 @@ }, { "cell_type": "markdown", - "id": "1080bf06", + "id": "eeefaf91", "metadata": { "editable": true }, @@ -2183,7 +2183,7 @@ }, { "cell_type": "markdown", - "id": "030ff44d", + "id": "5573b029", "metadata": { "editable": true }, @@ -2198,7 +2198,7 @@ }, { "cell_type": "markdown", - "id": "69f5722e", + "id": "4b23842f", "metadata": { "editable": true }, @@ -2212,7 +2212,7 @@ }, { "cell_type": "markdown", - "id": "1c05f1a1", + "id": "dd0076b0", "metadata": { "editable": true }, @@ -2223,7 +2223,7 @@ { "cell_type": "code", "execution_count": 9, - "id": "7f2203f8", + "id": "7efda7c7", "metadata": { "collapsed": false, "editable": true @@ -2284,7 +2284,7 @@ }, { "cell_type": "markdown", - "id": "e20918f8", + "id": "3fbe8c78", "metadata": { "editable": true }, @@ -2306,7 +2306,7 @@ }, { "cell_type": "markdown", - "id": "c4cd2dc7", + "id": "c19e4fc9", "metadata": { "editable": true }, @@ -2319,7 +2319,7 @@ { "cell_type": "code", "execution_count": 10, - "id": "4d62141d", + "id": "a6152435", "metadata": { "collapsed": false, "editable": true @@ -2385,7 +2385,7 @@ }, { "cell_type": "markdown", - "id": "f6af1edb", + "id": "e09ac4c8", "metadata": { "editable": true }, @@ -2396,7 +2396,7 @@ { "cell_type": "code", "execution_count": 11, - "id": "69121b13", + "id": "bbe6496f", "metadata": { "collapsed": false, "editable": true @@ -2470,7 +2470,7 @@ }, { "cell_type": "markdown", - "id": "f7ace6c9", + "id": "c5da561f", "metadata": { "editable": true }, @@ -2482,7 +2482,7 @@ }, { "cell_type": "markdown", - "id": "9ec31a1d", + "id": "8a534cf4", "metadata": { "editable": true }, @@ -2503,7 +2503,7 @@ }, { "cell_type": "markdown", - "id": "397dc9e2", + "id": "92d7a080", "metadata": { "editable": true }, @@ -2535,7 +2535,7 @@ }, { "cell_type": "markdown", - "id": "9e34ba26", + "id": "ee681442", "metadata": { "editable": true }, @@ -2552,7 +2552,7 @@ }, { "cell_type": "markdown", - "id": "3192866d", + "id": "72c1fa11", "metadata": { "editable": true }, @@ -2565,7 +2565,7 @@ }, { "cell_type": "markdown", - "id": "d0eacb6e", + "id": "71d38dd0", "metadata": { "editable": true }, @@ -2578,7 +2578,7 @@ }, { "cell_type": "markdown", - "id": "7f71e422", + "id": "b7bb3ead", "metadata": { "editable": true }, @@ -2591,7 +2591,7 @@ }, { "cell_type": "markdown", - "id": "68fe139a", + "id": "97d8d740", "metadata": { "editable": true }, @@ -2605,7 +2605,7 @@ }, { "cell_type": "markdown", - "id": "1b27dbae", + "id": "bc0338a7", "metadata": { "editable": true }, @@ -2627,7 +2627,7 @@ }, { "cell_type": "markdown", - "id": "d3aabc29", + "id": "37958c74", "metadata": { "editable": true }, @@ -2642,7 +2642,7 @@ }, { "cell_type": "markdown", - "id": "a393e4dc", + "id": "3ccbe5c6", "metadata": { "editable": true }, @@ -2654,7 +2654,7 @@ }, { "cell_type": "markdown", - "id": "454663c9", + "id": "30f2e7f5", "metadata": { "editable": true }, @@ -2667,7 +2667,7 @@ }, { "cell_type": "markdown", - "id": "2a2a9b25", + "id": "7907627a", "metadata": { "editable": true }, @@ -2681,7 +2681,7 @@ }, { "cell_type": "markdown", - "id": "dd7bc0c2", + "id": "1548f978", "metadata": { "editable": true }, @@ -2692,7 +2692,7 @@ { "cell_type": "code", "execution_count": 12, - "id": "7fc19ed3", + "id": "1ff042d7", "metadata": { "collapsed": false, "editable": true @@ -2717,7 +2717,7 @@ }, { "cell_type": "markdown", - "id": "b33dfafd", + "id": "6432d4b2", "metadata": { "editable": true }, @@ -2733,7 +2733,7 @@ }, { "cell_type": "markdown", - "id": "224cb6c6", + "id": "10ab1703", "metadata": { "editable": true }, @@ -2754,7 +2754,7 @@ }, { "cell_type": "markdown", - "id": "eea96e98", + "id": "dcd27ba0", "metadata": { "editable": true }, @@ -2774,7 +2774,7 @@ }, { "cell_type": "markdown", - "id": "91dcd74e", + "id": "7821fa5c", "metadata": { "editable": true }, @@ -2793,7 +2793,7 @@ { "cell_type": "code", "execution_count": 13, - "id": "0de8404e", + "id": "614e2a05", "metadata": { "collapsed": false, "editable": true @@ -2828,7 +2828,7 @@ }, { "cell_type": "markdown", - "id": "17e6deb7", + "id": "a6ddafd4", "metadata": { "editable": true }, @@ -2841,7 +2841,7 @@ { "cell_type": "code", "execution_count": 14, - "id": "21917f1d", + "id": "7dcb30ce", "metadata": { "collapsed": false, "editable": true @@ -2918,7 +2918,7 @@ }, { "cell_type": "markdown", - "id": "18e21a9b", + "id": "b6da28d8", "metadata": { "editable": true }, @@ -2933,7 +2933,7 @@ }, { "cell_type": "markdown", - "id": "90ea9013", + "id": "78015cb6", "metadata": { "editable": true }, @@ -2948,7 +2948,7 @@ }, { "cell_type": "markdown", - "id": "f9fa1bee", + "id": "ffc08e68", "metadata": { "editable": true }, @@ -2960,7 +2960,7 @@ }, { "cell_type": "markdown", - "id": "66018d26", + "id": "87a46d34", "metadata": { "editable": true }, @@ -2978,7 +2978,7 @@ }, { "cell_type": "markdown", - "id": "4e034652", + "id": "bc9a6e62", "metadata": { "editable": true }, @@ -2997,7 +2997,7 @@ }, { "cell_type": "markdown", - "id": "41863ff4", + "id": "4b91b7f6", "metadata": { "editable": true }, @@ -3009,7 +3009,7 @@ }, { "cell_type": "markdown", - "id": "dd86f43f", + "id": "7d7bad0b", "metadata": { "editable": true }, @@ -3019,7 +3019,7 @@ }, { "cell_type": "markdown", - "id": "be1bec74", + "id": "af03f313", "metadata": { "editable": true }, @@ -3035,7 +3035,7 @@ }, { "cell_type": "markdown", - "id": "f67f13cf", + "id": "3217e5f7", "metadata": { "editable": true }, @@ -3047,7 +3047,7 @@ }, { "cell_type": "markdown", - "id": "b00e4d9b", + "id": "fb12bd4a", "metadata": { "editable": true }, @@ -3057,7 +3057,7 @@ }, { "cell_type": "markdown", - "id": "645ea43b", + "id": "8d417735", "metadata": { "editable": true }, @@ -3069,7 +3069,7 @@ }, { "cell_type": "markdown", - "id": "f3155674", + "id": "7d6b3895", "metadata": { "editable": true }, @@ -3079,7 +3079,7 @@ }, { "cell_type": "markdown", - "id": "38c39b06", + "id": "691a6f0f", "metadata": { "editable": true }, @@ -3091,7 +3091,7 @@ }, { "cell_type": "markdown", - "id": "d09fe7eb", + "id": "6a5b7108", "metadata": { "editable": true }, @@ -3107,7 +3107,7 @@ }, { "cell_type": "markdown", - "id": "66fd53f5", + "id": "3cd05aba", "metadata": { "editable": true }, @@ -3119,7 +3119,7 @@ }, { "cell_type": "markdown", - "id": "1b552b2f", + "id": "d4792e11", "metadata": { "editable": true }, @@ -3152,7 +3152,7 @@ }, { "cell_type": "markdown", - "id": "bec97cf5", + "id": "a4814aed", "metadata": { "editable": true }, @@ -3164,7 +3164,7 @@ }, { "cell_type": "markdown", - "id": "d482e2fe", + "id": "9f289fcc", "metadata": { "editable": true }, @@ -3182,7 +3182,7 @@ }, { "cell_type": "markdown", - "id": "d0b6d3dc", + "id": "78417908", "metadata": { "editable": true }, @@ -3192,7 +3192,7 @@ }, { "cell_type": "markdown", - "id": "0510a924", + "id": "876f9f5e", "metadata": { "editable": true }, @@ -3223,7 +3223,7 @@ }, { "cell_type": "markdown", - "id": "05cebe24", + "id": "453c1007", "metadata": { "editable": true }, @@ -3238,7 +3238,7 @@ }, { "cell_type": "markdown", - "id": "f6feb0dd", + "id": "4fe0ca27", "metadata": { "editable": true }, @@ -3256,7 +3256,7 @@ }, { "cell_type": "markdown", - "id": "5c3d2f41", + "id": "b951e0f4", "metadata": { "editable": true }, @@ -3268,7 +3268,7 @@ }, { "cell_type": "markdown", - "id": "718af0d0", + "id": "b54a1702", "metadata": { "editable": true }, @@ -3280,7 +3280,7 @@ }, { "cell_type": "markdown", - "id": "34265c75", + "id": "24ec5c11", "metadata": { "editable": true }, @@ -3298,7 +3298,7 @@ }, { "cell_type": "markdown", - "id": "20ed6d82", + "id": "78f1652f", "metadata": { "editable": true }, @@ -3321,7 +3321,7 @@ }, { "cell_type": "markdown", - "id": "89caf637", + "id": "52633dca", "metadata": { "editable": true }, @@ -3339,7 +3339,7 @@ }, { "cell_type": "markdown", - "id": "f5e7a1ac", + "id": "37968950", "metadata": { "editable": true }, @@ -3351,7 +3351,7 @@ }, { "cell_type": "markdown", - "id": "5337c5ea", + "id": "51ab2eb9", "metadata": { "editable": true }, @@ -3363,7 +3363,7 @@ }, { "cell_type": "markdown", - "id": "9fa7d7c4", + "id": "eebe09f7", "metadata": { "editable": true }, @@ -3375,7 +3375,7 @@ }, { "cell_type": "markdown", - "id": "475f8da0", + "id": "0cccbdd6", "metadata": { "editable": true }, @@ -3387,7 +3387,7 @@ }, { "cell_type": "markdown", - "id": "8b337b55", + "id": "9662f3cc", "metadata": { "editable": true }, @@ -3399,7 +3399,7 @@ }, { "cell_type": "markdown", - "id": "71e429cc", + "id": "1aa8631f", "metadata": { "editable": true }, @@ -3416,7 +3416,7 @@ }, { "cell_type": "markdown", - "id": "391659f3", + "id": "e521cfb6", "metadata": { "editable": true }, @@ -3435,7 +3435,7 @@ }, { "cell_type": "markdown", - "id": "e1d8467c", + "id": "20fb3297", "metadata": { "editable": true }, @@ -3447,7 +3447,7 @@ }, { "cell_type": "markdown", - "id": "75e2c931", + "id": "a182d872", "metadata": { "editable": true }, @@ -3467,7 +3467,7 @@ }, { "cell_type": "markdown", - "id": "09303519", + "id": "a58a5b9c", "metadata": { "editable": true }, @@ -3505,7 +3505,7 @@ }, { "cell_type": "markdown", - "id": "f1be6f6b", + "id": "4fb66616", "metadata": { "editable": true }, @@ -3517,7 +3517,7 @@ }, { "cell_type": "markdown", - "id": "992f4a02", + "id": "a82f419d", "metadata": { "editable": true }, @@ -3527,7 +3527,7 @@ }, { "cell_type": "markdown", - "id": "a4b6d08e", + "id": "223e7cf3", "metadata": { "editable": true }, @@ -3539,7 +3539,7 @@ }, { "cell_type": "markdown", - "id": "c7d831ef", + "id": "871dd59d", "metadata": { "editable": true }, @@ -3550,7 +3550,7 @@ { "cell_type": "code", "execution_count": 15, - "id": "54acdf56", + "id": "0a5f00f5", "metadata": { "collapsed": false, "editable": true @@ -3595,7 +3595,7 @@ }, { "cell_type": "markdown", - "id": "9a753f05", + "id": "327a951f", "metadata": { "editable": true }, @@ -3612,7 +3612,7 @@ { "cell_type": "code", "execution_count": 16, - "id": "b8d502a7", + "id": "a7fb499f", "metadata": { "collapsed": false, "editable": true @@ -3640,7 +3640,7 @@ }, { "cell_type": "markdown", - "id": "db184594", + "id": "f3c96428", "metadata": { "editable": true }, @@ -3655,7 +3655,7 @@ { "cell_type": "code", "execution_count": 17, - "id": "61613ed6", + "id": "c9b33438", "metadata": { "collapsed": false, "editable": true @@ -3699,7 +3699,7 @@ }, { "cell_type": "markdown", - "id": "7543cea7", + "id": "564db542", "metadata": { "editable": true }, @@ -3709,7 +3709,7 @@ }, { "cell_type": "markdown", - "id": "a6c094d3", + "id": "4438a690", "metadata": { "editable": true }, @@ -3720,7 +3720,7 @@ { "cell_type": "code", "execution_count": 18, - "id": "e57abf0c", + "id": "d23e4aba", "metadata": { "collapsed": false, "editable": true @@ -3748,7 +3748,7 @@ }, { "cell_type": "markdown", - "id": "f63bb9f5", + "id": "f8ab8d5a", "metadata": { "editable": true }, @@ -3763,7 +3763,7 @@ }, { "cell_type": "markdown", - "id": "11623240", + "id": "baa23ffa", "metadata": { "editable": true }, @@ -3774,7 +3774,7 @@ { "cell_type": "code", "execution_count": 19, - "id": "148dce0d", + "id": "319def2d", "metadata": { "collapsed": false, "editable": true @@ -3802,7 +3802,7 @@ }, { "cell_type": "markdown", - "id": "7890ae18", + "id": "e1dd86c6", "metadata": { "editable": true }, @@ -3813,7 +3813,7 @@ { "cell_type": "code", "execution_count": 20, - "id": "1ec6199c", + "id": "e757e084", "metadata": { "collapsed": false, "editable": true @@ -3838,7 +3838,7 @@ }, { "cell_type": "markdown", - "id": "e87650cc", + "id": "eebfe3dd", "metadata": { "editable": true }, @@ -3849,7 +3849,7 @@ { "cell_type": "code", "execution_count": 21, - "id": "57766f6a", + "id": "c4830d3a", "metadata": { "collapsed": false, "editable": true @@ -3885,7 +3885,7 @@ { "cell_type": "code", "execution_count": 22, - "id": "43f2b332", + "id": "dc55ee7b", "metadata": { "collapsed": false, "editable": true @@ -3905,7 +3905,7 @@ }, { "cell_type": "markdown", - "id": "b3defe0d", + "id": "6cb061b0", "metadata": { "editable": true }, @@ -3916,7 +3916,7 @@ { "cell_type": "code", "execution_count": 23, - "id": "666d028d", + "id": "e468cdc4", "metadata": { "collapsed": false, "editable": true @@ -3954,7 +3954,7 @@ }, { "cell_type": "markdown", - "id": "1e69976c", + "id": "679da7ef", "metadata": { "editable": true }, @@ -3964,7 +3964,7 @@ }, { "cell_type": "markdown", - "id": "f8ca253c", + "id": "0ce4bd6b", "metadata": { "editable": true }, @@ -3978,7 +3978,7 @@ { "cell_type": "code", "execution_count": 24, - "id": "6f717821", + "id": "9838c269", "metadata": { "collapsed": false, "editable": true @@ -4000,7 +4000,7 @@ }, { "cell_type": "markdown", - "id": "06c934c9", + "id": "4457eba3", "metadata": { "editable": true }, @@ -4010,7 +4010,7 @@ }, { "cell_type": "markdown", - "id": "d199bb84", + "id": "d6303490", "metadata": { "editable": true }, @@ -4021,7 +4021,7 @@ { "cell_type": "code", "execution_count": 25, - "id": "49073461", + "id": "7dd580f0", "metadata": { "collapsed": false, "editable": true @@ -4043,7 +4043,7 @@ }, { "cell_type": "markdown", - "id": "673157aa", + "id": "27629408", "metadata": { "editable": true }, @@ -4056,7 +4056,7 @@ { "cell_type": "code", "execution_count": 26, - "id": "3aecc3d4", + "id": "30c2ae6c", "metadata": { "collapsed": false, "editable": true @@ -4081,7 +4081,7 @@ }, { "cell_type": "markdown", - "id": "5eab280c", + "id": "106e6e20", "metadata": { "editable": true }, @@ -4093,7 +4093,7 @@ { "cell_type": "code", "execution_count": 27, - "id": "ddff4eb6", + "id": "4aa28359", "metadata": { "collapsed": false, "editable": true @@ -4108,7 +4108,7 @@ }, { "cell_type": "markdown", - "id": "ebfad0f4", + "id": "0b41db61", "metadata": { "editable": true }, @@ -4123,7 +4123,7 @@ { "cell_type": "code", "execution_count": 28, - "id": "d99a8175", + "id": "0a3e77af", "metadata": { "collapsed": false, "editable": true @@ -4183,7 +4183,7 @@ }, { "cell_type": "markdown", - "id": "36efc0e8", + "id": "d66b1248", "metadata": { "editable": true }, @@ -4194,7 +4194,7 @@ { "cell_type": "code", "execution_count": 29, - "id": "7f807aeb", + "id": "62fd741e", "metadata": { "collapsed": false, "editable": true @@ -4258,7 +4258,67 @@ }, { "cell_type": "markdown", - "id": "02751092", + "id": "ea73ead4", + "metadata": { + "editable": true + }, + "source": [ + "## But noen of these can compete with Newton's method" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "85769e6a", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "# Using Newton's method\n", + "from random import random, seed\n", + "import numpy as np\n", + "import autograd.numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from autograd import grad\n", + "\n", + "def CostOLS(beta):\n", + " return (1.0/n)*np.sum((y-X @ beta)**2)\n", + "\n", + "n = 100\n", + "x = 2*np.random.rand(n,1)\n", + "y = 4+3*x+np.random.randn(n,1)\n", + "\n", + "X = np.c_[np.ones((n,1)), x]\n", + "XT_X = X.T @ X\n", + "beta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)\n", + "print(\"Own inversion\")\n", + "print(beta_linreg)\n", + "# Hessian matrix\n", + "H = (2.0/n)* XT_X\n", + "# Note that here the Hessian does not depend on the parameters beta\n", + "invH = np.linalg.pinv(H)\n", + "EigValues, EigVectors = np.linalg.eig(H)\n", + "print(f\"Eigenvalues of Hessian Matrix:{EigValues}\")\n", + "\n", + "beta = np.random.randn(2,1)\n", + "Niterations = 5\n", + "\n", + "# define the gradient\n", + "training_gradient = grad(CostOLS)\n", + "\n", + "for iter in range(Niterations):\n", + " gradients = training_gradient(beta)\n", + " beta -= invH @ gradients\n", + " print(iter,gradients[0],gradients[1])\n", + "print(\"beta from own Newton code\")\n", + "print(beta)" + ] + }, + { + "cell_type": "markdown", + "id": "c72ceec5", "metadata": { "editable": true }, @@ -4269,8 +4329,8 @@ }, { "cell_type": "code", - "execution_count": 30, - "id": "c958f749", + "execution_count": 31, + "id": "03605486", "metadata": { "collapsed": false, "editable": true @@ -4354,7 +4414,7 @@ }, { "cell_type": "markdown", - "id": "177b5e84", + "id": "44a2e49e", "metadata": { "editable": true }, @@ -4364,8 +4424,8 @@ }, { "cell_type": "code", - "execution_count": 31, - "id": "755d4596", + "execution_count": 32, + "id": "605503fb", "metadata": { "collapsed": false, "editable": true @@ -4443,7 +4503,90 @@ }, { "cell_type": "markdown", - "id": "66103f4e", + "id": "423148e4", + "metadata": { + "editable": true + }, + "source": [ + "## Same problem but now with AdaGrad" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "b025c34c", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "# Using Autograd to calculate gradients using SGD\n", + "# OLS example\n", + "from random import random, seed\n", + "import numpy as np\n", + "import autograd.numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from autograd import grad\n", + "\n", + "# Note change from previous example\n", + "def CostOLS(y,X,theta):\n", + " return np.sum((y-X @ theta)**2)\n", + "\n", + "n = 100\n", + "x = 2*np.random.rand(n,1)\n", + "y = 4+3*x+np.random.randn(n,1)\n", + "\n", + "X = np.c_[np.ones((n,1)), x]\n", + "XT_X = X.T @ X\n", + "theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)\n", + "print(\"Own inversion\")\n", + "print(theta_linreg)\n", + "# Hessian matrix\n", + "H = (2.0/n)* XT_X\n", + "EigValues, EigVectors = np.linalg.eig(H)\n", + "print(f\"Eigenvalues of Hessian Matrix:{EigValues}\")\n", + "\n", + "theta = np.random.randn(2,1)\n", + "eta = 1.0/np.max(EigValues)\n", + "Niterations = 100\n", + "\n", + "# Note that we request the derivative wrt third argument (theta, 2 here)\n", + "training_gradient = grad(CostOLS,2)\n", + "\n", + "for iter in range(Niterations):\n", + " gradients = (1.0/n)*training_gradient(y, X, theta)\n", + " theta -= eta*gradients\n", + "print(\"theta from own gd\")\n", + "print(theta)\n", + "print(np.size(gradients))\n", + "\n", + "n_epochs = 50\n", + "M = 5 #size of each minibatch\n", + "m = int(n/M) #number of minibatches\n", + "theta = np.random.randn(2,1)\n", + "\n", + "\n", + "# Including AdaGrad\n", + "delta = 0.000001\n", + "for epoch in range(n_epochs):\n", + " for i in range(m):\n", + " random_index = M*np.random.randint(m)\n", + " xi = X[random_index:random_index+M]\n", + " yi = y[random_index:random_index+M]\n", + " gradients = (1.0/M)*training_gradient(yi, xi, theta)\n", + " # calculate squared gradient by Hadamard multiplication\n", + " r -= gradients*gradients\n", + " # compute update\n", + " update = (1.0/delta+np.sqrt(r))*gradients\n", + " theta = eta*update\n", + "print(\"theta from own AdaGrad\")\n", + "print(theta)" + ] + }, + { + "cell_type": "markdown", + "id": "a1aebdad", "metadata": { "editable": true }, @@ -4453,8 +4596,8 @@ }, { "cell_type": "code", - "execution_count": 32, - "id": "b1f5ba7c", + "execution_count": 34, + "id": "1091d4ae", "metadata": { "collapsed": false, "editable": true @@ -4498,7 +4641,7 @@ }, { "cell_type": "markdown", - "id": "152c46ee", + "id": "f08388a7", "metadata": { "editable": true }, @@ -4516,8 +4659,8 @@ }, { "cell_type": "code", - "execution_count": 33, - "id": "34a7c3a4", + "execution_count": 35, + "id": "131b959a", "metadata": { "collapsed": false, "editable": true @@ -4534,6 +4677,22 @@ "derivative_fn = grad(sum_logistic)\n", "print(derivative_fn(x_small))" ] + }, + { + "cell_type": "markdown", + "id": "426864d6", + "metadata": { + "editable": true + }, + "source": [ + "## Weekend challenge\n", + "\n", + "* Try to run the above codes and implement the stochastic gradient descent with RMSprop and ADAM.\n", + "\n", + "* Add a more complicated function and study the rate of convergence for the derivatives as function of the different methods\n", + "\n", + "* Extend from linear regression to logistic regression." + ] } ], "metadata": {}, diff --git a/doc/src/week39/adagrad.py b/doc/src/week39/adagrad.py new file mode 100644 index 000000000..783124245 --- /dev/null +++ b/doc/src/week39/adagrad.py @@ -0,0 +1,67 @@ +# 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 + +# Note change from previous example +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 = 100 + +# 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) +print(np.size(gradients)) + +n_epochs = 50 +M = 5 #size of each minibatch +m = int(n/M) #number of minibatches +theta = np.random.randn(2,1) + + +# Including AdaGrad +delta = 0.000001 +for epoch in range(n_epochs): + for i in range(m): + random_index = M*np.random.randint(m) + xi = X[random_index:random_index+M] + yi = y[random_index:random_index+M] + gradients = (1.0/M)*training_gradient(yi, xi, theta) + # calculate squared gradient by Hadamard multiplication + r -= gradients*gradients + # compute update + update = (1.0/delta+np.sqrt(r))*gradients + theta = eta*update +print("theta from own AdaGrad") +print(theta) + + + + + + diff --git a/doc/src/week39/test.py b/doc/src/week39/codes/momentumGD.py similarity index 100% rename from doc/src/week39/test.py rename to doc/src/week39/codes/momentumGD.py diff --git a/doc/src/week39/sgd.py b/doc/src/week39/codes/momentumSGD.py similarity index 100% rename from doc/src/week39/sgd.py rename to doc/src/week39/codes/momentumSGD.py diff --git a/doc/src/week39/codes/newton.py b/doc/src/week39/codes/newton.py new file mode 100644 index 000000000..3c1cb5f0f --- /dev/null +++ b/doc/src/week39/codes/newton.py @@ -0,0 +1,38 @@ +# Using Newton's method +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 +beta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(beta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +# Note that here the Hessian does not depend on the parameters beta +invH = np.linalg.pinv(H) +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +beta = np.random.randn(2,1) +Niterations = 5 + +# define the gradient +training_gradient = grad(CostOLS) + +for iter in range(Niterations): + gradients = training_gradient(beta) + beta -= invH @ gradients + print(iter,gradients[0],gradients[1]) +print("beta from own Newton code") +print(beta) diff --git a/doc/src/week39/test.do.txt b/doc/src/week39/test.do.txt deleted file mode 100644 index 52a49347a..000000000 --- a/doc/src/week39/test.do.txt +++ /dev/null @@ -1,293 +0,0 @@ -TITLE: Codes -AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University -DATE: today -!bc pycod -from numpy import asarray -from numpy import arange -from numpy.random import rand -from numpy.random import seed -from matplotlib import pyplot - -# objective function -def objective(x): - return x**2.0 - -# derivative of objective function -def derivative(x): - return x * 2.0 - -# gradient descent algorithm -def gradient_descent(objective, derivative, bounds, n_iter, step_size): - # track all solutions - solutions, scores = list(), list() - # generate an initial point - solution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0]) - # run the gradient descent - for i in range(n_iter): - # calculate gradient - gradient = derivative(solution) - # take a step - solution = solution - step_size * gradient - # evaluate candidate point - solution_eval = objective(solution) - # store solution - solutions.append(solution) - scores.append(solution_eval) - # report progress - print('>%d f(%s) = %.5f' % (i, solution, solution_eval)) - return [solutions, scores] - -# seed the pseudo random number generator -seed(4) -# define range for input -bounds = asarray([[-1.0, 1.0]]) -# define the total iterations -n_iter = 30 -# define the step size -step_size = 0.1 -# perform the gradient descent search -solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size) -# sample input range uniformly at 0.1 increments -inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1) -# compute targets -results = objective(inputs) -# create a line plot of input vs result -pyplot.plot(inputs, results) -# plot the solutions found -pyplot.plot(solutions, scores, '.-', color='red') -# show the plot -pyplot.show() - -!ec - - -!split -===== Same code but now with momentum gradient descent ===== - -!bc pycod -from numpy import asarray -from numpy import arange -from numpy.random import rand -from numpy.random import seed -from matplotlib import pyplot - -# objective function -def objective(x): - return x**2.0 - -# derivative of objective function -def derivative(x): - return x * 2.0 - -# gradient descent algorithm -def gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum): - # track all solutions - solutions, scores = list(), list() - # generate an initial point - solution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0]) - # keep track of the change - change = 0.0 - # run the gradient descent - for i in range(n_iter): - # calculate gradient - gradient = derivative(solution) - # calculate update - new_change = step_size * gradient + momentum * change - # take a step - solution = solution - new_change - # save the change - change = new_change - # evaluate candidate point - solution_eval = objective(solution) - # store solution - solutions.append(solution) - scores.append(solution_eval) - # report progress - print('>%d f(%s) = %.5f' % (i, solution, solution_eval)) - return [solutions, scores] - -# seed the pseudo random number generator -seed(4) -# define range for input -bounds = asarray([[-1.0, 1.0]]) -# define the total iterations -n_iter = 30 -# define the step size -step_size = 0.1 -# define momentum -momentum = 0.3 -# perform the gradient descent search with momentum -solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum) -# sample input range uniformly at 0.1 increments -inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1) -# compute targets -results = objective(inputs) -# create a line plot of input vs result -pyplot.plot(inputs, results) -# plot the solutions found -pyplot.plot(solutions, scores, '.-', color='red') -# show the plot -pyplot.show() -!ec - - - -!split -===== Using Autograd with OLS ===== - -We conclude the part on optmization by showing how we can make codes -for linear regression and logistic regression using _autograd_. The -first example shows results with ordinary leats squares. - -!bc pycod -# Using Autograd to calculate gradients for OLS -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 -# define the gradient -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() - -!ec - - -!split -===== Including Stochastic Gradient Descent with Autograd ===== -In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using _autograd_. - -!bc pycod -# 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 - -# Note change from previous example -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 = 100 - -# 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 = M*np.random.randint(m) - xi = X[random_index:random_index+M] - yi = y[random_index:random_index+M] - gradients = (1.0/M)*training_gradient(yi, xi, theta) - eta = learning_schedule(epoch*m+i) - theta = theta - eta*gradients -print("theta from own sdg") -print(theta) - - -!ec - - - - - -!bc pycod -import jax.numpy as jnp -from jax import grad, jit, vmap - -def sum_logistic(x): - return jnp.sum(1.0 / (1.0 + jnp.exp(-x))) - -x_small = jnp.arange(3.) -derivative_fn = grad(sum_logistic) -print(derivative_fn(x_small)) - -!ec diff --git a/doc/src/week39/week39.do.txt b/doc/src/week39/week39.do.txt index 705aae706..26a1b4311 100644 --- a/doc/src/week39/week39.do.txt +++ b/doc/src/week39/week39.do.txt @@ -2228,11 +2228,53 @@ for iter in range(Niterations): print("theta from own gd wth momentum") print(theta) - - - !ec +!split +===== But noen of these can compete with Newton's method ===== + +!bc pycod +# Using Newton's method +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 +beta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(beta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +# Note that here the Hessian does not depend on the parameters beta +invH = np.linalg.pinv(H) +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +beta = np.random.randn(2,1) +Niterations = 5 + +# define the gradient +training_gradient = grad(CostOLS) + +for iter in range(Niterations): + gradients = training_gradient(beta) + beta -= invH @ gradients + print(iter,gradients[0],gradients[1]) +print("beta from own Newton code") +print(beta) +!ec + + !split ===== Including Stochastic Gradient Descent with Autograd ===== In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using _autograd_. @@ -2389,11 +2431,72 @@ print(theta) !ec +!split +===== Same problem but now with AdaGrad ===== +!bc pycod +# 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 + +# Note change from previous example +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 = 100 + +# 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) +print(np.size(gradients)) + +n_epochs = 50 +M = 5 #size of each minibatch +m = int(n/M) #number of minibatches +theta = np.random.randn(2,1) +# Including AdaGrad +delta = 0.000001 +for epoch in range(n_epochs): + for i in range(m): + random_index = M*np.random.randint(m) + xi = X[random_index:random_index+M] + yi = y[random_index:random_index+M] + gradients = (1.0/M)*training_gradient(yi, xi, theta) + # calculate squared gradient by Hadamard multiplication + r -= gradients*gradients + # compute update + update = (1.0/delta+np.sqrt(r))*gradients + theta = eta*update +print("theta from own AdaGrad") +print(theta) - - +!ec !split @@ -2459,3 +2562,10 @@ derivative_fn = grad(sum_logistic) print(derivative_fn(x_small)) !ec + +!split +===== Weekend challenge ===== + +* Try to run the above codes and implement the stochastic gradient descent with RMSprop and ADAM. +* Add a more complicated function and study the rate of convergence for the derivatives as function of the different methods +* Extend from linear regression to logistic regression.