From 321d23bb80f0b8f5f0f0b624b8ec78143d1c7c8e Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Mon, 29 Sep 2025 07:14:51 +0200 Subject: [PATCH] Update week40.do.txt --- doc/src/week40/week40.do.txt | 1581 ++++++++++++++++------------------ 1 file changed, 734 insertions(+), 847 deletions(-) diff --git a/doc/src/week40/week40.do.txt b/doc/src/week40/week40.do.txt index f44cf251a..6770a0d48 100644 --- a/doc/src/week40/week40.do.txt +++ b/doc/src/week40/week40.do.txt @@ -8,7 +8,7 @@ DATE: September 29-October 3, 2025 ===== Lecture Monday September 30, 2024 ===== !bblock o Logistic regression and gradient descent, examples on how to code -o Automatic differentiation and gradient descent, examples using Logistic regression +#o Automatic differentiation and gradient descent, examples using Logistic regression o Start with the basics of Neural Networks, setting up the basic steps, from the simple perceptron model to the multi-layer perceptron model # o "Video of lecture":"https://youtu.be/jdJoOrCIdII" # o Whiteboard notes at URL:"https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2024/NotesSeptember30.pdf" @@ -17,13 +17,13 @@ o Start with the basics of Neural Networks, setting up the basic steps, from the !split ===== Suggested readings and videos ===== !bblock Readings and Videos: - o The lecture notes for week 40 (these notes) - o For a good discussion on gradient methods, we would like to recommend Goodfellow et al section 4.3-4.5 and sections 8.3-8.6. We will come back to the latter chapter in our discussion of Neural networks as well. - o For neural networks we recommend Goodfellow et al chapter 6 and Raschka et al chapter 2 (contains also material about gradient descent) and chapter 11 (we will use this next week) +o The lecture notes for week 40 (these notes) +# o For a good discussion on gradient methods, we would like to recommend Goodfellow et al section 4.3-4.5 and# sections 8.3-8.6. We will come back to the latter chapter in our discussion of Neural networks as well. +o For neural networks we recommend Goodfellow et al chapter 6 and Raschka et al chapter 2 (contains also material about gradient descent) and chapter 11 (we will use this next week) # o Video on gradient descent at URL:"https://www.youtube.com/watch?v=sDv4f4s2SB8" - o Video on automatic differentiation at URL:"https://www.youtube.com/watch?v=wG_nF1awSSY" - o Neural Networks demystified at URL:"https://www.youtube.com/watch?v=bxe2T-V8XRs&list=PLiaHhY2iBX9hdHaRr6b7XevZtgZRa1PoU&ab_channel=WelchLabs" - o Building Neural Networks from scratch at URL:https://www.youtube.com/watch?v=Wo5dMEP_BbI&list=PLQVvvaa0QuDcjD5BAw2DxE6OF2tius3V3&ab_channel=sentdex" +# o Video on automatic differentiation at URL:"https://www.youtube.com/watch?v=wG_nF1awSSY" +o Neural Networks demystified at URL:"https://www.youtube.com/watch?v=bxe2T-V8XRs&list=PLiaHhY2iBX9hdHaRr6b7XevZtgZRa1PoU&ab_channel=WelchLabs" +o Building Neural Networks from scratch at URL:https://www.youtube.com/watch?v=Wo5dMEP_BbI&list=PLQVvvaa0QuDcjD5BAw2DxE6OF2tius3V3&ab_channel=sentdex" !eblock !split @@ -37,895 +37,781 @@ o Start with the basics of Neural Networks, setting up the basic steps, from the -!split -===== Automatic differentiation ===== - -"Automatic differentiation (AD)":"https://en.wikipedia.org/wiki/Automatic_differentiation", -also called algorithmic -differentiation or computational differentiation,is a set of -techniques to numerically evaluate the derivative of a function -specified by a computer program. AD exploits the fact that every -computer program, no matter how complicated, executes a sequence of -elementary arithmetic operations (addition, subtraction, -multiplication, division, etc.) and elementary functions (exp, log, -sin, cos, etc.). By applying the chain rule repeatedly to these -operations, derivatives of arbitrary order can be computed -automatically, accurately to working precision, and using at most a -small constant factor more arithmetic operations than the original -program. - -Automatic differentiation is neither: - -* Symbolic differentiation, nor -* Numerical differentiation (the method of finite differences). - -Symbolic differentiation can lead to inefficient code and faces the -difficulty of converting a computer program into a single expression, -while numerical differentiation can introduce round-off errors in the -discretization process and cancellation - - - -Python has tools for so-called _automatic differentiation_. -Consider the following example -!bt -\[ -f(x) = \sin\left(2\pi x + x^2\right) -\] -!et -which has the following derivative -!bt -\[ -f'(x) = \cos\left(2\pi x + x^2\right)\left(2\pi + 2x\right) -\] -!et -Using _autograd_ we have - -!bc pycod -import autograd.numpy as np - -# To do elementwise differentiation: -from autograd import elementwise_grad as egrad - -# To plot: -import matplotlib.pyplot as plt - - -def f(x): - return np.sin(2*np.pi*x + x**2) - -def f_grad_analytic(x): - return np.cos(2*np.pi*x + x**2)*(2*np.pi + 2*x) - -# Do the comparison: -x = np.linspace(0,1,1000) - -f_grad = egrad(f) - -computed = f_grad(x) -analytic = f_grad_analytic(x) - -plt.title('Derivative computed from Autograd compared with the analytical derivative') -plt.plot(x,computed,label='autograd') -plt.plot(x,analytic,label='analytic') - -plt.xlabel('x') -plt.ylabel('y') -plt.legend() - -plt.show() - -print("The max absolute difference is: %g"%(np.max(np.abs(computed - analytic)))) -!ec !split -===== Using autograd ===== +===== Logistic Regression, from last week ===== -Here we -experiment with what kind of functions Autograd is capable -of finding the gradient of. The following Python functions are just -meant to illustrate what Autograd can do, but please feel free to -experiment with other, possibly more complicated, functions as well. - -!bc pycod -import autograd.numpy as np -from autograd import grad - -def f1(x): - return x**3 + 1 - -f1_grad = grad(f1) - -# Remember to send in float as argument to the computed gradient from Autograd! -a = 1.0 - -# See the evaluated gradient at a using autograd: -print("The gradient of f1 evaluated at a = %g using autograd is: %g"%(a,f1_grad(a))) - -# Compare with the analytical derivative, that is f1'(x) = 3*x**2 -grad_analytical = 3*a**2 -print("The gradient of f1 evaluated at a = %g by finding the analytic expression is: %g"%(a,grad_analytical)) -!ec - - -!split -===== Autograd with more complicated functions ===== - -To differentiate with respect to two (or more) arguments of a Python -function, Autograd need to know at which variable the function if -being differentiated with respect to. - -!bc pycod -import autograd.numpy as np -from autograd import grad -def f2(x1,x2): - return 3*x1**3 + x2*(x1 - 5) + 1 - -# By sending the argument 0, Autograd will compute the derivative w.r.t the first variable, in this case x1 -f2_grad_x1 = grad(f2,0) - -# ... and differentiate w.r.t x2 by sending 1 as an additional arugment to grad -f2_grad_x2 = grad(f2,1) - -x1 = 1.0 -x2 = 3.0 - -print("Evaluating at x1 = %g, x2 = %g"%(x1,x2)) -print("-"*30) - -# Compare with the analytical derivatives: - -# Derivative of f2 w.r.t x1 is: 9*x1**2 + x2: -f2_grad_x1_analytical = 9*x1**2 + x2 - -# Derivative of f2 w.r.t x2 is: x1 - 5: -f2_grad_x2_analytical = x1 - 5 - -# See the evaluated derivations: -print("The derivative of f2 w.r.t x1: %g"%( f2_grad_x1(x1,x2) )) -print("The analytical derivative of f2 w.r.t x1: %g"%( f2_grad_x1(x1,x2) )) - -print() - -print("The derivative of f2 w.r.t x2: %g"%( f2_grad_x2(x1,x2) )) -print("The analytical derivative of f2 w.r.t x2: %g"%( f2_grad_x2(x1,x2) )) -!ec - -Note that the grad function will not produce the true gradient of the function. The true gradient of a function with two or more variables will produce a vector, where each element is the function differentiated w.r.t a variable. - - -!split -===== More complicated functions using the elements of their arguments directly ===== - -!bc pycod -import autograd.numpy as np -from autograd import grad -def f3(x): # Assumes x is an array of length 5 or higher - return 2*x[0] + 3*x[1] + 5*x[2] + 7*x[3] + 11*x[4]**2 - -f3_grad = grad(f3) - -x = np.linspace(0,4,5) - -# Print the computed gradient: -print("The computed gradient of f3 is: ", f3_grad(x)) - -# The analytical gradient is: (2, 3, 5, 7, 22*x[4]) -f3_grad_analytical = np.array([2, 3, 5, 7, 22*x[4]]) - -# Print the analytical gradient: -print("The analytical gradient of f3 is: ", f3_grad_analytical) -!ec - -Note that in this case, when sending an array as input argument, the -output from Autograd is another array. This is the true gradient of -the function, as opposed to the function in the previous example. By -using arrays to represent the variables, the output from Autograd -might be easier to work with, as the output is closer to what one -could expect form a gradient-evaluting function. +In linear regression our main interest was centered on learning the +coefficients of a functional fit (say a polynomial) in order to be +able to predict the response of a continuous variable on some unseen +data. The fit to the continuous variable $y_i$ is based on some +independent variables $\bm{x}_i$. Linear regression resulted in +analytical expressions for standard ordinary Least Squares or Ridge +regression (in terms of matrices to invert) for several quantities, +ranging from the variance and thereby the confidence intervals of the +parameters $\bm{\theta}$ to the mean squared error. If we can invert +the product of the design matrices, linear regression gives then a +simple recipe for fitting our data. !split -===== Functions using mathematical functions from Numpy ===== +===== Classification problems ===== -!bc pycod -import autograd.numpy as np -from autograd import grad -def f4(x): - return np.sqrt(1+x**2) + np.exp(x) + np.sin(2*np.pi*x) -f4_grad = grad(f4) +Classification problems, however, are concerned with outcomes taking +the form of discrete variables (i.e. categories). We may for example, +on the basis of DNA sequencing for a number of patients, like to find +out which mutations are important for a certain disease; or based on +scans of various patients' brains, figure out if there is a tumor or +not; or given a specific physical system, we'd like to identify its +state, say whether it is an ordered or disordered system (typical +situation in solid state physics); or classify the status of a +patient, whether she/he has a stroke or not and many other similar +situations. -x = 2.7 +The most common situation we encounter when we apply logistic +regression is that of two possible outcomes, normally denoted as a +binary outcome, true or false, positive or negative, success or +failure etc. -# Print the computed derivative: -print("The computed derivative of f4 at x = %g is: %g"%(x,f4_grad(x))) +!split +===== Optimization and Deep learning ===== + +Logistic regression will also serve as our stepping stone towards +neural network algorithms and supervised deep learning. For logistic +learning, the minimization of the cost function leads to a non-linear +equation in the parameters $\bm{\theta}$. The optimization of the +problem calls therefore for minimization algorithms. + +As we have discussed earlier, this forms the +bottle neck of all machine learning algorithms, namely how to find +reliable minima of a multi-variable function. This leads us to the +family of gradient descent methods. The latter are the working horses +of basically all modern machine learning algorithms. + +We note also that many of the topics discussed here on logistic +regression are also commonly used in modern supervised Deep Learning +models, as we will see later. + + +!split +===== Basics ===== + +We consider the case where the outputs/targets, also called the +responses or the outcomes, $y_i$ are discrete and only take values +from $k=0,\dots,K-1$ (i.e. $K$ classes). + +The goal is to predict the +output classes from the design matrix $\bm{X}\in\mathbb{R}^{n\times p}$ +made of $n$ samples, each of which carries $p$ features or predictors. The +primary goal is to identify the classes to which new unseen samples +belong. + +Last week we specialized to the case of two classes only, with outputs +$y_i=0$ and $y_i=1$. Our outcomes could represent the status of a +credit card user that could default or not on her/his credit card +debt. That is + + +!bt +\[ +y_i = \begin{bmatrix} 0 & \mathrm{no}\\ 1 & \mathrm{yes} \end{bmatrix}. +\] +!et -# The analytical derivative is: x/sqrt(1 + x**2) + exp(x) + cos(2*pi*x)*2*pi -f4_grad_analytical = x/np.sqrt(1 + x**2) + np.exp(x) + np.cos(2*np.pi*x)*2*np.pi -# Print the analytical gradient: -print("The analytical gradient of f4 at x = %g is: %g"%(x,f4_grad_analytical)) -!ec !split -===== More autograd ===== +===== Two parameters ===== -!bc pycod -import autograd.numpy as np -from autograd import grad -def f5(x): - if x >= 0: - return x**2 - else: - return -3*x + 1 +We assume now that we have two classes with $y_i$ either $0$ or $1$. Furthermore we assume also that we have only two parameters $\theta$ in our fitting of the Sigmoid function, that is we define probabilities +!bt +\begin{align*} +p(y_i=1|x_i,\bm{\theta}) &= \frac{\exp{(\theta_0+\theta_1x_i)}}{1+\exp{(\theta_0+\theta_1x_i)}},\nonumber\\ +p(y_i=0|x_i,\bm{\theta}) &= 1 - p(y_i=1|x_i,\bm{\theta}), +\end{align*} +!et +where $\bm{\theta}$ are the weights we wish to extract from data, in our case $\theta_0$ and $\theta_1$. -f5_grad = grad(f5) +Note that we used +!bt +\[ +p(y_i=0\vert x_i, \bm{\theta}) = 1-p(y_i=1\vert x_i, \bm{\theta}). +\] +!et -x = 2.7 +!split +===== Maximum likelihood ===== -# Print the computed derivative: -print("The computed derivative of f5 at x = %g is: %g"%(x,f5_grad(x))) -!ec +In order to define the total likelihood for all possible outcomes from a +dataset $\mathcal{D}=\{(y_i,x_i)\}$, with the binary labels +$y_i\in\{0,1\}$ and where the data points are drawn independently, we use the so-called "Maximum Likelihood Estimation":"https://en.wikipedia.org/wiki/Maximum_likelihood_estimation" (MLE) principle. +We aim thus at maximizing +the probability of seeing the observed data. We can then approximate the +likelihood in terms of the product of the individual probabilities of a specific outcome $y_i$, that is +!bt +\begin{align*} +P(\mathcal{D}|\bm{\theta})& = \prod_{i=1}^n \left[p(y_i=1|x_i,\bm{\theta})\right]^{y_i}\left[1-p(y_i=1|x_i,\bm{\theta}))\right]^{1-y_i}\nonumber \\ +\end{align*} +!et +from which we obtain the log-likelihood and our _cost/loss_ function +!bt +\[ +\mathcal{C}(\bm{\theta}) = \sum_{i=1}^n \left( y_i\log{p(y_i=1|x_i,\bm{\theta})} + (1-y_i)\log\left[1-p(y_i=1|x_i,\bm{\theta}))\right]\right). +\] +!et + +!split +===== The cost function rewritten ===== + +Reordering the logarithms, we can rewrite the _cost/loss_ function as +!bt +\[ +\mathcal{C}(\bm{\theta}) = \sum_{i=1}^n \left(y_i(\theta_0+\theta_1x_i) -\log{(1+\exp{(\theta_0+\theta_1x_i)})}\right). +\] +!et + +The maximum likelihood estimator is defined as the set of parameters that maximize the log-likelihood where we maximize with respect to $\theta$. +Since the cost (error) function is just the negative log-likelihood, for logistic regression we have that +!bt +\[ +\mathcal{C}(\bm{\theta})=-\sum_{i=1}^n \left(y_i(\theta_0+\theta_1x_i) -\log{(1+\exp{(\theta_0+\theta_1x_i)})}\right). +\] +!et +This equation is known in statistics as the _cross entropy_. Finally, we note that just as in linear regression, +in practice we often supplement the cross-entropy with additional regularization terms, usually $L_1$ and $L_2$ regularization as we did for Ridge and Lasso regression. + +!split +===== Minimizing the cross entropy ===== + +The cross entropy is a convex function of the weights $\bm{\theta}$ and, +therefore, any local minimizer is a global minimizer. + + +Minimizing this +cost function with respect to the two parameters $\theta_0$ and $\theta_1$ we obtain + +!bt +\[ +\frac{\partial \mathcal{C}(\bm{\theta})}{\partial \theta_0} = -\sum_{i=1}^n \left(y_i -\frac{\exp{(\theta_0+\theta_1x_i)}}{1+\exp{(\theta_0+\theta_1x_i)}}\right), +\] +!et +and +!bt +\[ +\frac{\partial \mathcal{C}(\bm{\theta})}{\partial \theta_1} = -\sum_{i=1}^n \left(y_ix_i -x_i\frac{\exp{(\theta_0+\theta_1x_i)}}{1+\exp{(\theta_0+\theta_1x_i)}}\right). +\] +!et + +!split +===== A more compact expression ===== + +Let us now define a vector $\bm{y}$ with $n$ elements $y_i$, an +$n\times p$ matrix $\bm{X}$ which contains the $x_i$ values and a +vector $\bm{p}$ of fitted probabilities $p(y_i\vert x_i,\bm{\theta})$. We can rewrite in a more compact form the first +derivative of the cost function as + +!bt +\[ +\frac{\partial \mathcal{C}(\bm{\theta})}{\partial \bm{\theta}} = -\bm{X}^T\left(\bm{y}-\bm{p}\right). +\] +!et + +If we in addition define a diagonal matrix $\bm{W}$ with elements +$p(y_i\vert x_i,\bm{\theta})(1-p(y_i\vert x_i,\bm{\theta})$, we can obtain a compact expression of the second derivative as + +!bt +\[ +\frac{\partial^2 \mathcal{C}(\bm{\theta})}{\partial \bm{\theta}\partial \bm{\theta}^T} = \bm{X}^T\bm{W}\bm{X}. +\] +!et + +!split +===== Extending to more predictors ===== + +Within a binary classification problem, we can easily expand our model to include multiple predictors. Our ratio between likelihoods is then with $p$ predictors +!bt +\[ +\log{ \frac{p(\bm{\theta}\bm{x})}{1-p(\bm{\theta}\bm{x})}} = \theta_0+\theta_1x_1+\theta_2x_2+\dots+\theta_px_p. +\] +!et +Here we defined $\bm{x}=[1,x_1,x_2,\dots,x_p]$ and $\bm{\theta}=[\theta_0, \theta_1, \dots, \theta_p]$ leading to +!bt +\[ +p(\bm{\theta}\bm{x})=\frac{ \exp{(\theta_0+\theta_1x_1+\theta_2x_2+\dots+\theta_px_p)}}{1+\exp{(\theta_0+\theta_1x_1+\theta_2x_2+\dots+\theta_px_p)}}. +\] +!et + +!split +===== Including more classes ===== + +Till now we have mainly focused on two classes, the so-called binary +system. Suppose we wish to extend to $K$ classes. Let us for the sake +of simplicity assume we have only two predictors. We have then following model + +!bt +\[ +\log{\frac{p(C=1\vert x)}{p(K\vert x)}} = \theta_{10}+\theta_{11}x_1, +\] +!et +and +!bt +\[ +\log{\frac{p(C=2\vert x)}{p(K\vert x)}} = \theta_{20}+\theta_{21}x_1, +\] +!et +and so on till the class $C=K-1$ class +!bt +\[ +\log{\frac{p(C=K-1\vert x)}{p(K\vert x)}} = \theta_{(K-1)0}+\theta_{(K-1)1}x_1, +\] +!et + +and the model is specified in term of $K-1$ so-called log-odds or +_logit_ transformations. !split -===== And with loops ===== +===== More classes ===== -!bc pycod -import autograd.numpy as np -from autograd import grad -def f6_for(x): - val = 0 - for i in range(10): - val = val + x**i - return val +In our discussion of neural networks we will encounter the above again +in terms of a slightly modified function, the so-called _Softmax_ function. -def f6_while(x): - val = 0 - i = 0 - while i < 10: - val = val + x**i - i = i + 1 - return val +The softmax function is used in various multiclass classification +methods, such as multinomial logistic regression (also known as +softmax regression), multiclass linear discriminant analysis, naive +Bayes classifiers, and artificial neural networks. Specifically, in +multinomial logistic regression and linear discriminant analysis, the +input to the function is the result of $K$ distinct linear functions, +and the predicted probability for the $k$-th class given a sample +vector $\bm{x}$ and a weighting vector $\bm{\theta}$ is (with two +predictors): -f6_for_grad = grad(f6_for) -f6_while_grad = grad(f6_while) +!bt +\[ +p(C=k\vert \mathbf {x} )=\frac{\exp{(\theta_{k0}+\theta_{k1}x_1)}}{1+\sum_{l=1}^{K-1}\exp{(\theta_{l0}+\theta_{l1}x_1)}}. +\] +!et +It is easy to extend to more predictors. The final class is +!bt +\[ +p(C=K\vert \mathbf {x} )=\frac{1}{1+\sum_{l=1}^{K-1}\exp{(\theta_{l0}+\theta_{l1}x_1)}}, +\] +!et -x = 0.5 +and they sum to one. Our earlier discussions were all specialized to +the case with two classes only. It is easy to see from the above that +what we derived earlier is compatible with these equations. -# Print the computed derivaties of f6_for and f6_while -print("The computed derivative of f6_for at x = %g is: %g"%(x,f6_for_grad(x))) -print("The computed derivative of f6_while at x = %g is: %g"%(x,f6_while_grad(x))) -!ec -!bc pycod -import autograd.numpy as np -from autograd import grad -# Both of the functions are implementation of the sum: sum(x**i) for i = 0, ..., 9 -# The analytical derivative is: sum(i*x**(i-1)) -f6_grad_analytical = 0 -for i in range(10): - f6_grad_analytical += i*x**(i-1) +To find the optimal parameters we would typically use a gradient +descent method. Newton's method and gradient descent methods are +discussed in the material on "optimization +methods":"https://compphysics.github.io/MachineLearning/doc/pub/Splines/html/Splines-bs.html". -print("The analytical derivative of f6 at x = %g is: %g"%(x,f6_grad_analytical)) -!ec -!split -===== Using recursion ===== -!bc pycod -import autograd.numpy as np -from autograd import grad -def f7(n): # Assume that n is an integer - if n == 1 or n == 0: - return 1 - else: - return n*f7(n-1) - -f7_grad = grad(f7) - -n = 2.0 - -print("The computed derivative of f7 at n = %d is: %g"%(n,f7_grad(n))) - -# The function f7 is an implementation of the factorial of n. -# By using the product rule, one can find that the derivative is: - -f7_grad_analytical = 0 -for i in range(int(n)-1): - tmp = 1 - for k in range(int(n)-1): - if k != i: - tmp *= (n - k) - f7_grad_analytical += tmp - -print("The analytical derivative of f7 at n = %d is: %g"%(n,f7_grad_analytical)) - -!ec -Note that if n is equal to zero or one, Autograd will give an error message. This message appears when the output is independent on input. !split -===== Using Autograd with OLS ===== +===== Optimization, the central part of any Machine Learning algortithm ===== -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. +Almost every problem in machine learning and data science starts with +a dataset $X$, a model $g(\theta)$, which is a function of the +parameters $\theta$ and a cost function $C(X, g(\theta))$ that allows +us to judge how well the model $g(\theta)$ explains the observations +$X$. The model is fit by finding the values of $\theta$ that minimize +the cost function. Ideally we would be able to solve for $\theta$ +analytically, however this is not possible in general and we must use +some approximative/numerical method to compute the minimum. + + +!split +===== Revisiting our Logistic Regression case ===== + +In our discussion on Logistic Regression we studied the +case of +two classes, with $y_i$ either +$0$ or $1$. Furthermore we assumed also that we have only two +parameters $\theta$ in our fitting, that is we +defined probabilities + +!bt +\begin{align*} +p(y_i=1|x_i,\bm{\theta}) &= \frac{\exp{(\theta_0+\theta_1x_i)}}{1+\exp{(\theta_0+\theta_1x_i)}},\nonumber\\ +p(y_i=0|x_i,\bm{\theta}) &= 1 - p(y_i=1|x_i,\bm{\theta}), +\end{align*} +!et +where $\bm{\theta}$ are the weights we wish to extract from data, in our case $\theta_0$ and $\theta_1$. + +!split +===== The equations to solve ===== + +Our compact equations used a definition of a vector $\bm{y}$ with $n$ +elements $y_i$, an $n\times p$ matrix $\bm{X}$ which contains the +$x_i$ values and a vector $\bm{p}$ of fitted probabilities +$p(y_i\vert x_i,\bm{\theta})$. We rewrote in a more compact form +the first derivative of the cost function as + +!bt +\[ +\frac{\partial \mathcal{C}(\bm{\theta})}{\partial \bm{\theta}} = -\bm{X}^T\left(\bm{y}-\bm{p}\right). +\] +!et + +If we in addition define a diagonal matrix $\bm{W}$ with elements +$p(y_i\vert x_i,\bm{\theta})(1-p(y_i\vert x_i,\bm{\theta})$, we can obtain a compact expression of the second derivative as + +!bt +\[ +\frac{\partial^2 \mathcal{C}(\bm{\theta})}{\partial \bm{\theta}\partial \bm{\theta}^T} = \bm{X}^T\bm{W}\bm{X}. +\] +!et +This defines what is called the Hessian matrix. + +!split +===== Solving using Newton-Raphson's method ===== + +If we can set up these equations, Newton-Raphson's iterative method is normally the method of choice. It requires however that we can compute in an efficient way the matrices that define the first and second derivatives. + +Our iterative scheme is then given by + +!bt +\[ +\bm{\theta}^{\mathrm{new}} = \bm{\theta}^{\mathrm{old}}-\left(\frac{\partial^2 \mathcal{C}(\bm{\theta})}{\partial \bm{\theta}\partial \bm{\theta}^T}\right)^{-1}_{\bm{\theta}^{\mathrm{old}}}\times \left(\frac{\partial \mathcal{C}(\bm{\theta})}{\partial \bm{\theta}}\right)_{\bm{\theta}^{\mathrm{old}}}, +\] +!et +or in matrix form as + +!bt +\[ +\bm{\theta}^{\mathrm{new}} = \bm{\theta}^{\mathrm{old}}-\left(\bm{X}^T\bm{W}\bm{X} \right)^{-1}\times \left(-\bm{X}^T(\bm{y}-\bm{p}) \right)_{\bm{\theta}^{\mathrm{old}}}. +\] +!et +The right-hand side is computed with the old values of $\theta$. + +If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement. + + + +!split +===== Example code for Logistic Regression ===== + +Here we make a class for Logistic regression. The code uses a simple data set and includes both a binary case and a multiclass case. !bc pycod -# Using Autograd to calculate gradients for OLS -from random import random, seed import numpy as np -import autograd.numpy as np + +class LogisticRegression: + """ + Logistic Regression for binary and multiclass classification. + """ + def __init__(self, lr=0.01, epochs=1000, fit_intercept=True, verbose=False): + self.lr = lr # Learning rate for gradient descent + self.epochs = epochs # Number of iterations + self.fit_intercept = fit_intercept # Whether to add intercept (bias) + self.verbose = verbose # Print loss during training if True + self.weights = None + self.multi_class = False # Will be determined at fit time + + def _add_intercept(self, X): + """Add intercept term (column of ones) to feature matrix.""" + intercept = np.ones((X.shape[0], 1)) + return np.concatenate((intercept, X), axis=1) + + def _sigmoid(self, z): + """Sigmoid function for binary logistic.""" + return 1 / (1 + np.exp(-z)) + + def _softmax(self, Z): + """Softmax function for multiclass logistic.""" + exp_Z = np.exp(Z - np.max(Z, axis=1, keepdims=True)) + return exp_Z / np.sum(exp_Z, axis=1, keepdims=True) + + def fit(self, X, y): + """ + Train the logistic regression model using gradient descent. + Supports binary (sigmoid) and multiclass (softmax) based on y. + """ + X = np.array(X) + y = np.array(y) + n_samples, n_features = X.shape + + # Add intercept if needed + if self.fit_intercept: + X = self._add_intercept(X) + n_features += 1 + + # Determine classes and mode (binary vs multiclass) + unique_classes = np.unique(y) + if len(unique_classes) > 2: + self.multi_class = True + else: + self.multi_class = False + + # ----- Multiclass case ----- + if self.multi_class: + n_classes = len(unique_classes) + # Map original labels to 0...n_classes-1 + class_to_index = {c: idx for idx, c in enumerate(unique_classes)} + y_indices = np.array([class_to_index[c] for c in y]) + # Initialize weight matrix (features x classes) + self.weights = np.zeros((n_features, n_classes)) + + # One-hot encode y + Y_onehot = np.zeros((n_samples, n_classes)) + Y_onehot[np.arange(n_samples), y_indices] = 1 + + # Gradient descent + for epoch in range(self.epochs): + scores = X.dot(self.weights) # Linear scores (n_samples x n_classes) + probs = self._softmax(scores) # Probabilities (n_samples x n_classes) + # Compute gradient (features x classes) + gradient = (1 / n_samples) * X.T.dot(probs - Y_onehot) + # Update weights + self.weights -= self.lr * gradient + + if self.verbose and epoch % 100 == 0: + # Compute current loss (categorical cross-entropy) + loss = -np.sum(Y_onehot * np.log(probs + 1e-15)) / n_samples + print(f"[Epoch {epoch}] Multiclass loss: {loss:.4f}") + + # ----- Binary case ----- + else: + # Convert y to 0/1 if not already + if not np.array_equal(unique_classes, [0, 1]): + # Map the two classes to 0 and 1 + class0, class1 = unique_classes + y_binary = np.where(y == class1, 1, 0) + else: + y_binary = y.copy().astype(int) + + # Initialize weights vector (features,) + self.weights = np.zeros(n_features) + + # Gradient descent + for epoch in range(self.epochs): + linear_model = X.dot(self.weights) # (n_samples,) + probs = self._sigmoid(linear_model) # (n_samples,) + # Gradient for binary cross-entropy + gradient = (1 / n_samples) * X.T.dot(probs - y_binary) + self.weights -= self.lr * gradient + + if self.verbose and epoch % 100 == 0: + # Compute binary cross-entropy loss + loss = -np.mean( + y_binary * np.log(probs + 1e-15) + + (1 - y_binary) * np.log(1 - probs + 1e-15) + ) + print(f"[Epoch {epoch}] Binary loss: {loss:.4f}") + + def predict_prob(self, X): + """ + Compute probability estimates. Returns a 1D array for binary or + a 2D array (n_samples x n_classes) for multiclass. + """ + X = np.array(X) + # Add intercept if the model used it + if self.fit_intercept: + X = self._add_intercept(X) + scores = X.dot(self.weights) + if self.multi_class: + return self._softmax(scores) + else: + return self._sigmoid(scores) + + def predict(self, X): + """ + Predict class labels for samples in X. + Returns integer class labels (0,1 for binary, or 0...C-1 for multiclass). + """ + probs = self.predict_prob(X) + if self.multi_class: + # Choose class with highest probability + return np.argmax(probs, axis=1) + else: + # Threshold at 0.5 for binary + return (probs >= 0.5).astype(int) +!ec + + +The class implements the sigmoid and softmax internally. During fit(), +we check the number of classes: if more than 2, we set +self.multi_class=True and perform multinomial logistic regression. We +one-hot encode the target vector and update a weight matrix with +softmax probabilities. Otherwise, we do standard binary logistic +regression, converting labels to 0/1 if needed and updating a weight +vector. In both cases we use batch gradient descent on the +cross-entropy loss (we add a small epsilon 1e-15 to logs for numerical +stability). Progress (loss) can be printed if verbose=True. + +!bc pycod +# Evaluation Metrics +#We define helper functions for accuracy and cross-entropy loss. Accuracy is the fraction of correct predictions . For loss, we compute the appropriate cross-entropy: + +def accuracy_score(y_true, y_pred): + """Accuracy = (# correct predictions) / (total samples).""" + y_true = np.array(y_true) + y_pred = np.array(y_pred) + return np.mean(y_true == y_pred) + +def binary_cross_entropy(y_true, y_prob): + """ + Binary cross-entropy loss. + y_true: true binary labels (0 or 1), y_prob: predicted probabilities for class 1. + """ + y_true = np.array(y_true) + y_prob = np.clip(np.array(y_prob), 1e-15, 1-1e-15) + return -np.mean(y_true * np.log(y_prob) + (1 - y_true) * np.log(1 - y_prob)) + +def categorical_cross_entropy(y_true, y_prob): + """ + Categorical cross-entropy loss for multiclass. + y_true: true labels (0...C-1), y_prob: array of predicted probabilities (n_samples x C). + """ + y_true = np.array(y_true, dtype=int) + y_prob = np.clip(np.array(y_prob), 1e-15, 1-1e-15) + # One-hot encode true labels + n_samples, n_classes = y_prob.shape + one_hot = np.zeros_like(y_prob) + one_hot[np.arange(n_samples), y_true] = 1 + # Compute cross-entropy + loss_vec = -np.sum(one_hot * np.log(y_prob), axis=1) + return np.mean(loss_vec) +!ec + + +=== Synthetic data generation === + +Binary classification data: Create two Gaussian clusters in 2D. For example, class 0 around mean [-2,-2] and class 1 around [2,2]. +Multiclass data: Create several Gaussian clusters (one per class) spread out in feature space. + + +!bc pycod +import numpy as np + +def generate_binary_data(n_samples=100, n_features=2, random_state=None): + """ + Generate synthetic binary classification data. + Returns (X, y) where X is (n_samples x n_features), y in {0,1}. + """ + rng = np.random.RandomState(random_state) + # Half samples for class 0, half for class 1 + n0 = n_samples // 2 + n1 = n_samples - n0 + # Class 0 around mean -2, class 1 around +2 + mean0 = -2 * np.ones(n_features) + mean1 = 2 * np.ones(n_features) + X0 = rng.randn(n0, n_features) + mean0 + X1 = rng.randn(n1, n_features) + mean1 + X = np.vstack((X0, X1)) + y = np.array([0]*n0 + [1]*n1) + return X, y + +def generate_multiclass_data(n_samples=150, n_features=2, n_classes=3, random_state=None): + """ + Generate synthetic multiclass data with n_classes Gaussian clusters. + """ + rng = np.random.RandomState(random_state) + X = [] + y = [] + samples_per_class = n_samples // n_classes + for cls in range(n_classes): + # Random cluster center for each class + center = rng.uniform(-5, 5, size=n_features) + Xi = rng.randn(samples_per_class, n_features) + center + yi = [cls] * samples_per_class + X.append(Xi) + y.extend(yi) + X = np.vstack(X) + y = np.array(y) + return X, y + + +# Generate and test on binary data +X_bin, y_bin = generate_binary_data(n_samples=200, n_features=2, random_state=42) +model_bin = LogisticRegression(lr=0.1, epochs=1000) +model_bin.fit(X_bin, y_bin) +y_prob_bin = model_bin.predict_prob(X_bin) # probabilities for class 1 +y_pred_bin = model_bin.predict(X_bin) # predicted classes 0 or 1 + +acc_bin = accuracy_score(y_bin, y_pred_bin) +loss_bin = binary_cross_entropy(y_bin, y_prob_bin) +print(f"Binary Classification - Accuracy: {acc_bin:.2f}, Cross-Entropy Loss: {loss_bin:.2f}") +#For multiclass: +# Generate and test on multiclass data +X_multi, y_multi = generate_multiclass_data(n_samples=300, n_features=2, n_classes=3, random_state=1) +model_multi = LogisticRegression(lr=0.1, epochs=1000) +model_multi.fit(X_multi, y_multi) +y_prob_multi = model_multi.predict_prob(X_multi) # (n_samples x 3) probabilities +y_pred_multi = model_multi.predict(X_multi) # predicted labels 0,1,2 + +acc_multi = accuracy_score(y_multi, y_pred_multi) +loss_multi = categorical_cross_entropy(y_multi, y_prob_multi) +print(f"Multiclass Classification - Accuracy: {acc_multi:.2f}, Cross-Entropy Loss: {loss_multi:.2f}") + +# CSV Export +import csv + +# Export binary results +with open('binary_results.csv', mode='w', newline='') as f: + writer = csv.writer(f) + writer.writerow(["TrueLabel", "PredictedLabel"]) + for true, pred in zip(y_bin, y_pred_bin): + writer.writerow([true, pred]) + +# Export multiclass results +with open('multiclass_results.csv', mode='w', newline='') as f: + writer = csv.writer(f) + writer.writerow(["TrueLabel", "PredictedLabel"]) + for true, pred in zip(y_multi, y_pred_multi): + writer.writerow([true, pred]) + +!ec + + +!split +===== Using _Scikit-learn_ ===== + +We show here how we can use a logistic regression case on a data set +included in _scikit_learn_, the so-called Wisconsin breast cancer data +using Logistic regression as our algorithm for classification. This is +a widely studied data set and can easily be included in demonstrations +of classification problems. + + +!bc pycod import matplotlib.pyplot as plt -from autograd import grad +import numpy as np +from sklearn.model_selection import train_test_split +from sklearn.datasets import load_breast_cancer +from sklearn.linear_model import LogisticRegression -def CostOLS(beta): - return (1.0/n)*np.sum((y-X @ beta)**2) +# Load the data +cancer = load_breast_cancer() -n = 100 -x = 2*np.random.rand(n,1) -y = 4+3*x+np.random.randn(n,1) +X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0) +print(X_train.shape) +print(X_test.shape) +# Logistic Regression +logreg = LogisticRegression(solver='lbfgs') +logreg.fit(X_train, y_train) +print("Test set accuracy with Logistic Regression: {:.2f}".format(logreg.score(X_test,y_test))) +!ec -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}") +!split +===== Using the correlation matrix ===== -theta = np.random.randn(2,1) -eta = 1.0/np.max(EigValues) -Niterations = 1000 -# define the gradient -training_gradient = grad(CostOLS) +In addition to the above scores, we could also study the covariance (and the correlation matrix). +We use _Pandas_ to compute the correlation matrix. +!bc pycod +import matplotlib.pyplot as plt +import numpy as np +from sklearn.model_selection import train_test_split +from sklearn.datasets import load_breast_cancer +from sklearn.linear_model import LogisticRegression +cancer = load_breast_cancer() +import pandas as pd +# Making a data frame +cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names) -for iter in range(Niterations): - gradients = training_gradient(theta) - theta -= eta*gradients -print("theta from own gd") -print(theta) +fig, axes = plt.subplots(15,2,figsize=(10,20)) +malignant = cancer.data[cancer.target == 0] +benign = cancer.data[cancer.target == 1] +ax = axes.ravel() -xnew = np.array([[0],[2]]) -Xnew = np.c_[np.ones((2,1)), xnew] -ypredict = Xnew.dot(theta) -ypredict2 = Xnew.dot(theta_linreg) +for i in range(30): + _, bins = np.histogram(cancer.data[:,i], bins =50) + ax[i].hist(malignant[:,i], bins = bins, alpha = 0.5) + ax[i].hist(benign[:,i], bins = bins, alpha = 0.5) + ax[i].set_title(cancer.feature_names[i]) + ax[i].set_yticks(()) +ax[0].set_xlabel("Feature magnitude") +ax[0].set_ylabel("Frequency") +ax[0].legend(["Malignant", "Benign"], loc ="best") +fig.tight_layout() +plt.show() -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 ') +import seaborn as sns +correlation_matrix = cancerpd.corr().round(1) +# use the heatmap function from seaborn to plot the correlation matrix +# annot = True to print the values inside the square +plt.figure(figsize=(15,8)) +sns.heatmap(data=correlation_matrix, annot=True) +plt.show() + + +!ec + +!split +===== Discussing the correlation data ===== + +In the above example we note two things. In the first plot we display +the overlap of benign and malignant tumors as functions of the various +features in the Wisconsin data set. We see that for +some of the features we can distinguish clearly the benign and +malignant cases while for other features we cannot. This can point to +us which features may be of greater interest when we wish to classify +a benign or not benign tumour. + +In the second figure we have computed the so-called correlation +matrix, which in our case with thirty features becomes a $30\times 30$ +matrix. + +We constructed this matrix using _pandas_ via the statements +!bc pycod +cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names) +!ec +and then +!bc pycod +correlation_matrix = cancerpd.corr().round(1) +!ec + +Diagonalizing this matrix we can in turn say something about which +features are of relevance and which are not. This leads us to +the classical Principal Component Analysis (PCA) theorem with +applications. This will be discussed later this semester. + + + +!split +===== Other measures in classification studies ===== +!bc pycod +import matplotlib.pyplot as plt +import numpy as np +from sklearn.model_selection import train_test_split +from sklearn.datasets import load_breast_cancer +from sklearn.linear_model import LogisticRegression + +# Load the data +cancer = load_breast_cancer() + +X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0) +print(X_train.shape) +print(X_test.shape) +# Logistic Regression +logreg = LogisticRegression(solver='lbfgs') +logreg.fit(X_train, y_train) + +from sklearn.preprocessing import LabelEncoder +from sklearn.model_selection import cross_validate +#Cross validation +accuracy = cross_validate(logreg,X_test,y_test,cv=10)['test_score'] +print(accuracy) +print("Test set accuracy with Logistic Regression: {:.2f}".format(logreg.score(X_test,y_test))) + +import scikitplot as skplt +y_pred = logreg.predict(X_test) +skplt.metrics.plot_confusion_matrix(y_test, y_pred, normalize=True) +plt.show() +y_probas = logreg.predict_proba(X_test) +skplt.metrics.plot_roc(y_test, y_probas) +plt.show() +skplt.metrics.plot_cumulative_gain(y_test, y_probas) plt.show() !ec -!split -===== Same code but now with momentum gradient descent ===== -!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 = 30 - -# define the gradient -training_gradient = grad(CostOLS) - -for iter in range(Niterations): - gradients = training_gradient(theta) - theta -= eta*gradients - print(iter,gradients[0],gradients[1]) -print("theta from own gd") -print(theta) - -# Now improve with momentum gradient descent -change = 0.0 -delta_momentum = 0.3 -for iter in range(Niterations): - # calculate gradient - gradients = training_gradient(theta) - # calculate update - new_change = eta*gradients+delta_momentum*change - # take a step - theta -= new_change - # save the change - change = new_change - print(iter,gradients[0],gradients[1]) -print("theta from own gd wth momentum") -print(theta) - -!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 = 1000 - -# Note that we request the derivative wrt third argument (theta, 2 here) -training_gradient = grad(CostOLS,2) - -for iter in range(Niterations): - gradients = (1.0/n)*training_gradient(y, X, theta) - theta -= eta*gradients -print("theta from own gd") -print(theta) - -xnew = np.array([[0],[2]]) -Xnew = np.c_[np.ones((2,1)), xnew] -ypredict = Xnew.dot(theta) -ypredict2 = Xnew.dot(theta_linreg) - -plt.plot(xnew, ypredict, "r-") -plt.plot(xnew, ypredict2, "b-") -plt.plot(x, y ,'ro') -plt.axis([0,2.0,0, 15.0]) -plt.xlabel(r'$x$') -plt.ylabel(r'$y$') -plt.title(r'Random numbers ') -plt.show() - -n_epochs = 50 -M = 5 #size of each minibatch -m = int(n/M) #number of minibatches -t0, t1 = 5, 50 -def learning_schedule(t): - return t0/(t+t1) - -theta = np.random.randn(2,1) - -for epoch in range(n_epochs): -# Can you figure out a better way of setting up the contributions to each batch? - for i in range(m): - random_index = 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 - - -!split -===== Same code but now with momentum gradient descent ===== -!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) - - -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) - -change = 0.0 -delta_momentum = 0.3 - -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) - eta = learning_schedule(epoch*m+i) - # calculate update - new_change = eta*gradients+delta_momentum*change - # take a step - theta -= new_change - # save the change - change = new_change -print("theta from own sdg with momentum") -print(theta) -!ec - - -!split -===== Similar (second order function now) problem but now with AdaGrad ===== -!bc pycod -# Using Autograd to calculate gradients using AdaGrad and Stochastic Gradient descent -# 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 = 1000 -x = np.random.rand(n,1) -y = 2.0+3*x +4*x*x - -X = np.c_[np.ones((n,1)), x, x*x] -XT_X = X.T @ X -theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) -print("Own inversion") -print(theta_linreg) - - -# Note that we request the derivative wrt third argument (theta, 2 here) -training_gradient = grad(CostOLS,2) -# Define parameters for Stochastic Gradient Descent -n_epochs = 50 -M = 5 #size of each minibatch -m = int(n/M) #number of minibatches -# Guess for unknown parameters theta -theta = np.random.randn(3,1) - -# Value for learning rate -eta = 0.01 -# Including AdaGrad parameter to avoid possible division by zero -delta = 1e-8 -for epoch in range(n_epochs): - Giter = 0.0 - 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) - Giter += gradients*gradients - update = gradients*eta/(delta+np.sqrt(Giter)) - theta -= update -print("theta from own AdaGrad") -print(theta) - - -!ec - -Running this code we note an almost perfect agreement with the results from matrix inversion. - -!split -===== RMSprop for adaptive learning rate with Stochastic Gradient Descent ===== -!bc pycod -# Using Autograd to calculate gradients using RMSprop and Stochastic Gradient descent -# 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 = 1000 -x = np.random.rand(n,1) -y = 2.0+3*x +4*x*x# +np.random.randn(n,1) - -X = np.c_[np.ones((n,1)), x, x*x] -XT_X = X.T @ X -theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) -print("Own inversion") -print(theta_linreg) - - -# Note that we request the derivative wrt third argument (theta, 2 here) -training_gradient = grad(CostOLS,2) -# Define parameters for Stochastic Gradient Descent -n_epochs = 50 -M = 5 #size of each minibatch -m = int(n/M) #number of minibatches -# Guess for unknown parameters theta -theta = np.random.randn(3,1) - -# Value for learning rate -eta = 0.01 -# Value for parameter rho -rho = 0.99 -# Including AdaGrad parameter to avoid possible division by zero -delta = 1e-8 -for epoch in range(n_epochs): - Giter = 0.0 - 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) - # Accumulated gradient - # Scaling with rho the new and the previous results - Giter = (rho*Giter+(1-rho)*gradients*gradients) - # Taking the diagonal only and inverting - update = gradients*eta/(delta+np.sqrt(Giter)) - # Hadamard product - theta -= update -print("theta from own RMSprop") -print(theta) -!ec - -!split -===== And finally "ADAM":"https://arxiv.org/pdf/1412.6980.pdf" ===== - -!bc pycod -# Using Autograd to calculate gradients using RMSprop and Stochastic Gradient descent -# 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 = 1000 -x = np.random.rand(n,1) -y = 2.0+3*x +4*x*x# +np.random.randn(n,1) - -X = np.c_[np.ones((n,1)), x, x*x] -XT_X = X.T @ X -theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) -print("Own inversion") -print(theta_linreg) - - -# Note that we request the derivative wrt third argument (theta, 2 here) -training_gradient = grad(CostOLS,2) -# Define parameters for Stochastic Gradient Descent -n_epochs = 50 -M = 5 #size of each minibatch -m = int(n/M) #number of minibatches -# Guess for unknown parameters theta -theta = np.random.randn(3,1) - -# Value for learning rate -eta = 0.01 -# Value for parameters beta1 and beta2, see https://arxiv.org/abs/1412.6980 -beta1 = 0.9 -beta2 = 0.999 -# Including AdaGrad parameter to avoid possible division by zero -delta = 1e-7 -iter = 0 -for epoch in range(n_epochs): - first_moment = 0.0 - second_moment = 0.0 - iter += 1 - 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) - # Computing moments first - first_moment = beta1*first_moment + (1-beta1)*gradients - second_moment = beta2*second_moment+(1-beta2)*gradients*gradients - first_term = first_moment/(1.0-beta1**iter) - second_term = second_moment/(1.0-beta2**iter) - # Scaling with rho the new and the previous results - update = eta*first_term/(np.sqrt(second_term)+delta) - theta -= update -print("theta from own ADAM") -print(theta) -!ec - -!split -===== And Logistic Regression ===== - -!bc pycod -import autograd.numpy as np -from autograd import grad - -def sigmoid(x): - return 0.5 * (np.tanh(x / 2.) + 1) - -def logistic_predictions(weights, inputs): - # Outputs probability of a label being true according to logistic model. - return sigmoid(np.dot(inputs, weights)) - -def training_loss(weights): - # Training loss is the negative log-likelihood of the training labels. - preds = logistic_predictions(weights, inputs) - label_probabilities = preds * targets + (1 - preds) * (1 - targets) - return -np.sum(np.log(label_probabilities)) - -# Build a toy dataset. -inputs = np.array([[0.52, 1.12, 0.77], - [0.88, -1.08, 0.15], - [0.52, 0.06, -1.30], - [0.74, -2.49, 1.39]]) -targets = np.array([True, True, False, True]) - -# Define a function that returns gradients of training loss using Autograd. -training_gradient_fun = grad(training_loss) - -# Optimize weights using gradient descent. -weights = np.array([0.0, 0.0, 0.0]) -print("Initial loss:", training_loss(weights)) -for i in range(100): - weights -= training_gradient_fun(weights) * 0.01 - -print("Trained loss:", training_loss(weights)) -!ec - - - - -===== Introducing "JAX":"https://jax.readthedocs.io/en/latest/" ===== - -Presently, instead of using _autograd_, we recommend using "JAX":"https://jax.readthedocs.io/en/latest/" - -_JAX_ is Autograd and "XLA (Accelerated Linear Algebra))":"https://www.tensorflow.org/xla", -brought together for high-performance numerical computing and machine learning research. -It provides composable transformations of Python+NumPy programs: differentiate, vectorize, parallelize, Just-In-Time compile to GPU/TPU, and more. - -=== Getting started with Jax, note the way we import numpy === -!bc pycod -import jax -import jax.numpy as jnp -import numpy as np -import matplotlib.pyplot as plt - -from jax import grad as jax_grad -!ec - - -=== A warm-up example === - -!bc pycod -def function(x): - return x**2 - -def analytical_gradient(x): - return 2*x - -def gradient_descent(starting_point, learning_rate, num_iterations, solver="analytical"): - x = starting_point - trajectory_x = [x] - trajectory_y = [function(x)] - - if solver == "analytical": - grad = analytical_gradient - elif solver == "jax": - grad = jax_grad(function) - x = jnp.float64(x) - learning_rate = jnp.float64(learning_rate) - - for _ in range(num_iterations): - - x = x - learning_rate * grad(x) - trajectory_x.append(x) - trajectory_y.append(function(x)) - - return trajectory_x, trajectory_y - -x = np.linspace(-5, 5, 100) -plt.plot(x, function(x), label="f(x)") - -descent_x, descent_y = gradient_descent(5, 0.1, 10, solver="analytical") -jax_descend_x, jax_descend_y = gradient_descent(5, 0.1, 10, solver="jax") - -plt.plot(descent_x, descent_y, label="Gradient descent", marker="o") -plt.plot(jax_descend_x, jax_descend_y, label="JAX", marker="x") -!ec - -=== A more advanced example === - -!bc pycod -backend = np - -def function(x): - return x*backend.sin(x**2 + 1) - -def analytical_gradient(x): - return backend.sin(x**2 + 1) + 2*x**2*backend.cos(x**2 + 1) - - -x = np.linspace(-5, 5, 100) -plt.plot(x, function(x), label="f(x)") - -descent_x, descent_y = gradient_descent(1, 0.01, 300, solver="analytical") - -# Change the backend to JAX -backend = jnp -jax_descend_x, jax_descend_y = gradient_descent(1, 0.01, 300, solver="jax") - -plt.scatter(descent_x, descent_y, label="Gradient descent", marker="v", s=10, color="red") -plt.scatter(jax_descend_x, jax_descend_y, label="JAX", marker="x", s=5, color="black") -!ec @@ -1542,3 +1428,4 @@ plt.show() !ec +