diff --git a/doc/src/week41/figures/simplenn1.png b/doc/src/week41/figures/simplenn1.png
new file mode 100644
index 000000000..3c87aa3ee
Binary files /dev/null and b/doc/src/week41/figures/simplenn1.png differ
diff --git a/doc/src/week41/figures/simplenn2.png b/doc/src/week41/figures/simplenn2.png
new file mode 100644
index 000000000..2ce83dd53
Binary files /dev/null and b/doc/src/week41/figures/simplenn2.png differ
diff --git a/doc/src/week41/week3.do.txt b/doc/src/week41/week3.do.txt
new file mode 100644
index 000000000..9480a011f
--- /dev/null
+++ b/doc/src/week41/week3.do.txt
@@ -0,0 +1,896 @@
+TITLE: January 29-February 2 : Advanced machine learning and data analysis for the physical sciences
+AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics and Center for Computing in Science Education, University of Oslo, Norway & Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University, East Lansing, Michigan, USA
+DATE: January 30
+
+
+
+!split
+===== Overview of third  week  =====
+
+!bblock 
+  o Discussion of possible projects
+  o Review of neural networks and automatic differentiation
+  o Discussion of codes
+  o "Video of lecture":"https://youtu.be/OUTFo0oJadU"
+  o "Link to material for project suggestions":"https://github.com/CompPhysics/AdvancedMachineLearning/tree/main/doc/Projects/2024/ProjectProposals"
+
+!eblock
+
+
+!split
+===== Mathematics of deep learning =====
+
+!bblock Two recent books online
+o "The Modern Mathematics of Deep Learning, by Julius Berner, Philipp Grohs, Gitta Kutyniok, Philipp Petersen":"https://arxiv.org/abs/2105.04026", published as "Mathematical Aspects of Deep Learning, pp. 1-111. Cambridge University Press, 2022":"https://doi.org/10.1017/9781009025096.002"
+
+o "Mathematical Introduction to Deep Learning: Methods, Implementations, and Theory, Arnulf Jentzen, Benno Kuckuck, Philippe von Wurstemberger":"https://doi.org/10.48550/arXiv.2310.20360"
+!eblock
+
+
+!split
+===== Reminder on books with hands-on material and codes =====
+!bblock
+* "Sebastian Rashcka et al, Machine learning with Sickit-Learn and PyTorch":"https://sebastianraschka.com/blog/2022/ml-pytorch-book.html"
+* "David Foster, Generative Deep Learning with TensorFlow":"https://www.oreilly.com/library/view/generative-deep-learning/9781098134174/ch01.html"
+* "Bali and Gavras, Generative AI with Python and TensorFlow 2":"https://github.com/PacktPublishing/Hands-On-Generative-AI-with-Python-and-TensorFlow-2"
+!eblock
+
+All three books have GitHub addresses from where  one can download all codes. We will borrow most of the material from these three texts as well as 
+from Goodfellow, Bengio and Courville's text "Deep Learning":"https://www.deeplearningbook.org/"
+
+
+!split
+===== Reading recommendations =====
+
+o Rashkca et al., chapter 11, jupyter-notebook sent separately, from "GitHub":"https://github.com/rasbt/machine-learning-book"
+o Goodfellow et al, chapter 6 and 7 contain most of the neural network background.
+
+!split
+===== Mathematics of deep learning and neural networks =====
+
+
+Neural networks, in its so-called feed-forward form, where each
+iterations contains a feed-forward stage and a back-propgagation
+stage, consist of series of affine matrix-matrix and matrix-vector
+multiplications. The unknown parameters (the so-called biases and
+weights which deternine the architecture of a neural network), are
+uptaded iteratively using the so-called back-propagation algorithm.
+This algorithm corresponds to the so-called reverse mode of 
+automatic differentation. 
+
+!split
+===== Basics of an NN =====
+
+A neural network consists of a series of hidden layers, in addition to
+the input and output layers.  Each layer $l$ has a set of parameters
+$\bm{\Theta}^{(l)}=(\bm{W}^{(l)},\bm{b}^{(l)})$ which are related to the
+parameters in other layers through a series of affine transformations,
+for a standard NN these are matrix-matrix and matrix-vector
+multiplications.  For all layers we will simply use a collective variable $\bm{\Theta}$.
+
+It consist of two basic steps:
+o a feed forward stage which takes a given input and produces a final output which is compared with the target values through our cost/loss function.
+o a back-propagation state where the unknown parameters $\bm{\Theta}$ are updated through the optimization of the their gradients. The expressions for the gradients are obtained via the chain rule, starting from the derivative of the cost/function.
+
+These two steps make up one iteration. This iterative process is continued till we reach an eventual stopping criterion.
+
+
+!split
+===== Overarching view of a neural network =====
+
+The architecture of a neural network defines our model. This model
+aims at describing some function $f(\bm{x}$ which represents
+some final result (outputs or tagrget values) given a specific inpput
+$\bm{x}$. Note that here $\bm{y}$ and $\bm{x}$ are not limited to be
+vectors.
+
+The architecture consists of
+o An input and an output layer where the input layer is defined by the inputs $\bm{x}$. The output layer produces the model ouput $\bm{\tilde{y}}$ which is compared with the target value $\bm{y}$
+o A given number of hidden layers and neurons/nodes/units for each layer (this may vary)
+o A given activation function $\sigma(\bm{z})$ with arguments $\bm{z}$ to be defined below. The activation functions may differ from layer to layer.
+o The last layer, normally called _output_ layer has normally an activation function tailored to the specific problem
+o Finally we define a so-called cost or loss function which is used to gauge the quality of our model. 
+
+
+!split
+===== The optimization problem =====
+
+The cost function is a function of the unknown parameters
+$\bm{\Theta}$ where the latter is a container for all possible
+parameters needed to define a neural network
+
+If we are dealing with a regression task a typical cost/loss function
+is the mean squared error
+!bt
+\[
+C(\bm{\Theta})=\frac{1}{n}\left\{\left(\bm{y}-\bm{X}\bm{\theta}\right)^T\left(\bm{y}-\bm{X}\bm{\theta}\right)\right\}.
+\]
+!et
+This function represents one of many possible ways to define
+the so-called cost function. Note that here we have assumed a linear dependence in terms of the paramters $\bm{\Theta}$. This is in general not the case.
+
+
+!split
+===== Parameters of neural networks =====
+For neural networks the parameters
+$\bm{\Theta}$ are given by the so-called weights and biases (to be
+defined below).
+
+The weights are given by matrix elements $w_{ij}^{(l)}$ where the
+superscript indicates the layer number. The biases are typically given
+by vector elements representing each single node of a given layer,
+that is $b_j^{(l)}$.
+
+!split
+===== Other ingredients of a neural network =====
+
+Having defined the architecture of a neural network, the optimization
+of the cost function with respect to the parameters $\bm{\Theta}$,
+involves the calculations of gradients and their optimization. The
+gradients represent the derivatives of a multidimensional object and
+are often approximated by various gradient methods, including
+o various quasi-Newton methods,
+o plain gradient descent (GD) with a constant learning rate $\eta$,
+o GD with momentum and other approximations to the learning rates such as
+  * Adapative gradient (ADAgrad)
+  * Root mean-square propagation (RMSprop)
+  * Adaptive gradient with momentum (ADAM) and many other
+o Stochastic gradient descent and various families of learning rate approximations
+
+!split
+===== Other parameters =====
+
+In addition to the above, there are often additional hyperparamaters
+which are included in the setup of a neural network. These will be
+discussed below.
+
+
+!split
+===== Universal approximation theorem =====
+
+The universal approximation theorem plays a central role in deep
+learning.  "Cybenko (1989)":"https://link.springer.com/article/10.1007/BF02551274" showed
+the following:
+
+!bblock
+Let $\sigma$ be any continuous sigmoidal function such that
+!bt
+\[
+\sigma(z) = \left\{\begin{array}{cc} 1 & z\rightarrow \infty\\ 0 & z \rightarrow -\infty \end{array}\right.
+\]
+!et
+Given a continuous and deterministic function $F(\bm{x})$ on the unit
+cube in $d$-dimensions $F\in [0,1]^d$, $x\in [0,1]^d$ and a parameter
+$\epsilon >0$, there is a one-layer (hidden) neural network
+$f(\bm{x};\bm{\Theta})$ with $\bm{\Theta}=(\bm{W},\bm{b})$ and $\bm{W}\in
+\mathbb{R}^{m\times n}$ and $\bm{b}\in \mathbb{R}^{n}$, for which
+!bt
+\[
+\vert F(\bm{x})-f(\bm{x};\bm{\Theta})\vert < \epsilon \hspace{0.1cm} \forall \bm{x}\in[0,1]^d.
+\]
+!et
+
+!eblock
+
+!split
+===== Some parallels from real analysis =====
+
+For those of you familiar with for example the "Stone-Weierstrass
+theorem":"https://en.wikipedia.org/wiki/Stone%E2%80%93Weierstrass_theorem"
+for polynomial approximations or the convergence criterion for Fourier
+series, there are similarities in the derivation of the proof for
+neural networks. 
+
+!split
+===== The approximation theorem in words =====
+
+_Any continuous function $y=F(\bm{x})$ supported on the unit cube in
+$d$-dimensions can be approximated by a one-layer sigmoidal network to
+arbitrary accuracy._
+
+"Hornik (1991)":"https://www.sciencedirect.com/science/article/abs/pii/089360809190009T" extended the theorem by letting any non-constant, bounded activation function to be included using that the expectation value
+!bt
+\[
+\mathbb{E}[\vert F(\bm{x})\vert^2] =\int_{\bm{x}\in D} \vert F(\bm{x})\vert^2p(\bm{x})d\bm{x} < \infty.
+\]
+!et
+Then we have
+!bt
+\[
+\mathbb{E}[\vert F(\bm{x})-f(\bm{x};\bm{\Theta})\vert^2] =\int_{\bm{x}\in D} \vert F(\bm{x})-f(\bm{x};\bm{\Theta})\vert^2p(\bm{x})d\bm{x} < \epsilon.
+\]
+!et
+
+
+!split
+===== More on the general approximation theorem =====
+
+None of the proofs give any insight into the relation between the
+number of of hidden layers and nodes and the approximation error
+$\epsilon$, nor the magnitudes of $\bm{W}$ and $\bm{b}$.
+
+Neural networks (NNs) have what we may call a kind of universality no matter what function we want to compute.
+
+!bblock
+It does not mean that an NN can be used to exactly compute any function. Rather, we get an approximation that is as good as we want. 
+!eblock
+
+!split
+===== Class of functions we can approximate =====
+
+!bblock
+The class of functions that can be approximated are the continuous ones.
+If the function $F(\bm{x})$ is discontinuous, it won't in general be possible to approximate it. However, an NN may still give an approximation even if we fail in some points.
+!eblock
+
+
+
+!split
+===== Setting up the equations for a neural network =====
+
+The questions we want to ask are how do changes in the biases and the
+weights in our network change the cost function and how can we use the
+final output to modify the weights and biases?
+
+To derive these equations let us start with a plain regression problem
+and define our cost function as
+
+!bt
+\[
+{\cal C}(\bm{\Theta})  =  \frac{1}{2}\sum_{i=1}^n\left(y_i - \tilde{y}_i\right)^2, 
+\]
+!et
+
+where the $y_i$s are our $n$ targets (the values we want to
+reproduce), while the outputs of the network after having propagated
+all inputs $\bm{x}$ are given by $\bm{\tilde{y}}_i$.
+
+
+!split
+===== Layout of a neural network with three hidden layers =====
+
+FIGURE: [figures/nn1.pdf, width=900 frac=1.0]
+
+!split
+===== Definitions =====
+
+With our definition of the targets $\bm{y}$, the outputs of the
+network $\bm{\tilde{y}}$ and the inputs $\bm{x}$ we
+define now the activation $z_j^l$ of node/neuron/unit $j$ of the
+$l$-th layer as a function of the bias, the weights which add up from
+the previous layer $l-1$ and the forward passes/outputs
+$\hat{a}^{l-1}$ from the previous layer as
+
+
+!bt
+\[
+z_j^l = \sum_{i=1}^{M_{l-1}}w_{ij}^la_i^{l-1}+b_j^l,
+\]
+!et
+
+where $b_k^l$ are the biases from layer $l$.  Here $M_{l-1}$
+represents the total number of nodes/neurons/units of layer $l-1$. The
+figure in the whiteboard notes illustrates this equation.  We can rewrite this in a more
+compact form as the matrix-vector products we discussed earlier,
+
+!bt
+\[
+\hat{z}^l = \left(\hat{W}^l\right)^T\hat{a}^{l-1}+\hat{b}^l.
+\]
+!et
+
+!split
+===== Inputs to the activation function =====
+
+With the activation values $\bm{z}^l$ we can in turn define the
+output of layer $l$ as $\bm{a}^l = f(\bm{z}^l)$ where $f$ is our
+activation function. In the examples here we will use the sigmoid
+function discussed in our logistic regression lectures. We will also use the same activation function $f$ for all layers
+and their nodes.  It means we have
+
+!bt
+\[
+a_j^l = \sigma(z_j^l) = \frac{1}{1+\exp{-(z_j^l)}}.
+\]
+!et
+
+
+
+
+!split
+===== Derivatives and the chain rule =====
+
+From the definition of the activation $z_j^l$ we have
+!bt
+\[
+\frac{\partial z_j^l}{\partial w_{ij}^l} = a_i^{l-1},
+\]
+!et
+and
+!bt
+\[
+\frac{\partial z_j^l}{\partial a_i^{l-1}} = w_{ji}^l. 
+\]
+!et
+
+With our definition of the activation function we have that (note that this function depends only on $z_j^l$)
+!bt
+\[
+\frac{\partial a_j^l}{\partial z_j^{l}} = a_j^l(1-a_j^l)=\sigma(z_j^l)(1-\sigma(z_j^l)). 
+\]
+!et
+
+
+!split
+===== Derivative of the cost function =====
+
+With these definitions we can now compute the derivative of the cost function in terms of the weights.
+
+Let us specialize to the output layer $l=L$. Our cost function is
+!bt
+\[
+{\cal C}(\bm{\Theta}^L)  =  \frac{1}{2}\sum_{i=1}^n\left(y_i - \tilde{y}_i\right)^2=\frac{1}{2}\sum_{i=1}^n\left(a_i^L - y_i\right)^2, 
+\]
+!et
+The derivative of this function with respect to the weights is
+
+!bt
+\[
+\frac{\partial{\cal C}(\bm{\Theta}^L)}{\partial w_{jk}^L}  =  \left(a_j^L - y_j\right)\frac{\partial a_j^L}{\partial w_{jk}^{L}}, 
+\]
+!et
+The last partial derivative can easily be computed and reads (by applying the chain rule)
+!bt
+\[
+\frac{\partial a_j^L}{\partial w_{jk}^{L}} = \frac{\partial a_j^L}{\partial z_{j}^{L}}\frac{\partial z_j^L}{\partial w_{jk}^{L}}=a_j^L(1-a_j^L)a_k^{L-1}.  
+\]
+!et
+
+
+
+!split
+===== Simpler examples first, and automatic differentiation =====
+
+In order to understand the back propagation algorithm and its
+derivation (an implementation of the chain rule), let us first digress
+with some simple examples. These examples are also meant to motivate
+the link with back propagation and "automatic differentiation":"https://en.wikipedia.org/wiki/Automatic_differentiation".
+
+!split
+===== Reminder on the chain rule and gradients =====
+
+If we have a multivariate function $f(x,y)$ where $x=x(t)$ and $y=y(t)$ are functions of a variable $t$, we have that the gradient of $f$ with respect to $t$ (without the explicit unit vector components)
+!bt
+\[
+\frac{df}{dt} = \begin{bmatrix}\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \end{bmatrix} \begin{bmatrix}\frac{\partial x}{\partial t} \\ \frac{\partial y}{\partial t} \end{bmatrix}=\frac{\partial f}{\partial x} \frac{\partial x}{\partial t} +\frac{\partial f}{\partial y} \frac{\partial y}{\partial t}. 
+\]
+!et
+
+
+!split
+===== Multivariable functions =====
+
+If we have a multivariate function $f(x,y)$ where $x=x(t,s)$ and $y=y(t,s)$ are functions of the variables $t$ and $s$, we have that the partial derivatives
+!bt
+\[
+\frac{\partial f}{\partial s}=\frac{\partial f}{\partial x}\frac{\partial x}{\partial s}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial s},
+\]
+!et
+and
+!bt
+\[
+\frac{\partial f}{\partial t}=\frac{\partial f}{\partial x}\frac{\partial x}{\partial t}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial t}.
+\]
+!et
+
+the gradient of $f$ with respect to $t$ and $s$ (without the explicit unit vector components)
+!bt
+\[
+\frac{df}{d(s,t)} = \begin{bmatrix}\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \end{bmatrix} \begin{bmatrix}\frac{\partial x}{\partial s}  &\frac{\partial x}{\partial t} \\ \frac{\partial y}{\partial s} & \frac{\partial y}{\partial t} \end{bmatrix}.
+\]
+!et
+
+
+
+!split
+===== Automatic differentiation through examples =====
+
+A great introduction to automatic differentiation is given by Baydin et al., see URL:"https://arxiv.org/abs/1502.05767".
+
+Automatic differentiation is a represented by a repeated application
+of the chain rule on well-known functions and allows for the
+calculation of derivatives to numerical precision. It is not the same
+as the calculation of symbolic derivatives via for example SymPy, nor
+does it use approximative formulae based on Taylor-expansions of a
+function around a given value. The latter are error prone due to
+truncation errors and values of the step size $\Delta$.
+
+!split
+===== Simple example =====
+
+Our first example is rather simple,
+!bt
+\[
+f(x) =\exp{x^2},
+\]
+!et
+with derivative
+!bt
+\[
+f'(x) =2x\exp{x^2}.
+\]
+!et
+We can use SymPy to extract the pertinent lines of Python code through the following simple example
+!bc pycod
+from __future__ import division
+from sympy import *
+x = symbols('x')
+expr = exp(x*x)
+simplify(expr)
+derivative = diff(expr,x)
+print(python(expr))
+print(python(derivative))
+!ec
+
+!split
+===== Smarter way of evaluating the above function =====
+If we study this function, we note that we can reduce the number of operations by introducing an intermediate variable
+!bt
+\[
+a = x^2,
+\]
+!et
+leading to 
+!bt
+\[
+f(x) = f(a(x)) = b= \exp{a}.
+\]
+!et
+
+We now assume that all operations can be counted in terms of equal
+floating point operations. This means that in order to calculate
+$f(x)$ we need first to square $x$ and then compute the exponential. We
+have thus two floating point operations only.
+
+!split
+===== Reducing the number of operations =====
+
+With the introduction of a precalculated quantity $a$ and thereby $f(x)$ we have that the derivative can be written as
+
+!bt
+\[
+f'(x) = 2xb,
+\]
+!et
+
+which reduces the number of operations from four in the orginal
+expression to two. This means that if we need to compute $f(x)$ and
+its derivative (a common task in optimizations), we have reduced the
+number of operations from six to four in total.
+
+_Note_ that the usage of a symbolic software like SymPy does not
+include such simplifications and the calculations of the function and
+the derivatives yield in general more floating point operations.
+
+!split
+===== Chain rule, forward and reverse modes =====
+
+In the above example we have introduced the variables $a$ and $b$, and our function is
+!bt
+\[
+f(x) = f(a(x)) = b= \exp{a},
+\]
+!et
+with $a=x^2$. We can decompose the derivative of $f$ with respect to $x$ as
+!bt
+\[
+\frac{df}{dx}=\frac{df}{db}\frac{db}{da}\frac{da}{dx}.
+\]
+!et
+
+We note that since $b=f(x)$ that
+!bt
+\[
+\frac{df}{db}=1,
+\]
+!et
+leading to
+!bt
+\[
+\frac{df}{dx}=\frac{db}{da}\frac{da}{dx}=2x\exp{x^2},
+\]
+!et
+as before.
+
+
+!split
+===== Forward and reverse modes =====
+
+We have that 
+!bt
+\[
+\frac{df}{dx}=\frac{df}{db}\frac{db}{da}\frac{da}{dx},
+\]
+!et
+which we can rewrite either as
+!bt
+\[
+\frac{df}{dx}=\left[\frac{df}{db}\frac{db}{da}\right]\frac{da}{dx},
+\]
+!et
+or
+!bt
+\[
+\frac{df}{dx}=\frac{df}{db}\left[\frac{db}{da}\frac{da}{dx}\right].
+\]
+!et
+
+The first expression is called reverse mode (or back propagation)
+since we start by evaluating the derivatives at the end point and then
+propagate backwards. This is the standard way of evaluating
+derivatives (gradients) when optimizing the parameters of a neural
+network.  In the context of deep learning this is computationally
+more efficient since the output of a neural network consists of either
+one or some few other output variables.
+
+The second equation defines the so-called  _forward mode_.
+
+
+!split
+===== More complicated function =====
+
+We increase our ambitions and introduce a slightly more complicated function
+!bt
+\[
+f(x) =\sqrt{x^2+exp{x^2}},
+\]
+!et
+with derivative
+!bt
+\[
+f'(x) =\frac{x(1+\exp{x^2})}{\sqrt{x^2+exp{x^2}}}.
+\]
+!et
+The corresponding SymPy code reads
+!bc pycod
+from __future__ import division
+from sympy import *
+x = symbols('x')
+expr = sqrt(x*x+exp(x*x))
+simplify(expr)
+derivative = diff(expr,x)
+print(python(expr))
+print(python(derivative))
+!ec
+
+!split
+===== Counting the number of floating point operations =====
+
+A simple count of operations shows that we need five operations for
+the function itself and ten for the derivative.  Fifteen operations in total if we wish to proceed with the above codes.
+
+Can we reduce this to
+say half the number of operations?
+
+!split
+=====  Defining intermediate operations =====
+
+We can indeed reduce the number of operation to half of those listed in the brute force approach above.
+We define the following quantities
+!bt
+\[
+a = x^2,
+\]
+!et
+and
+!bt
+\[
+b = \exp{x^2} = \exp{a},
+\]
+!et
+and
+!bt
+\[
+c= a+b,
+\]
+!et
+and
+!bt
+\[
+d=f(x)=\sqrt{c}.
+\]
+!et
+
+!split
+===== New expression for the derivative =====
+
+With these definitions we obtain the following partial derivatives 
+!bt
+\[
+\frac{\partial a}{\partial x} = 2x,
+\]
+!et
+and
+!bt
+\[
+\frac{\partial b}{\partial a} = \exp{a},
+\]
+!et
+and
+!bt
+\[
+\frac{\partial c}{\partial a} = 1,
+\]
+!et
+and
+!bt
+\[
+\frac{\partial c}{\partial b} = 1,
+\]
+!et
+and
+!bt
+\[
+\frac{\partial d}{\partial c} = \frac{1}{2\sqrt{c}},
+\]
+!et
+and finally
+!bt
+\[
+\frac{\partial f}{\partial d} = 1.
+\]
+!et
+
+!split
+===== Final derivatives =====
+Our final derivatives are thus
+!bt
+\[
+\frac{\partial f}{\partial c} = \frac{\partial f}{\partial d} \frac{\partial d}{\partial c}  = \frac{1}{2\sqrt{c}},
+\]
+!et
+!bt
+\[
+\frac{\partial f}{\partial b} = \frac{\partial f}{\partial c} \frac{\partial c}{\partial b}  = \frac{1}{2\sqrt{c}},
+\]
+!et
+!bt
+\[
+\frac{\partial f}{\partial a} = \frac{\partial f}{\partial c} \frac{\partial c}{\partial a}+
+\frac{\partial f}{\partial b} \frac{\partial b}{\partial a}  = \frac{1+\exp{a}}{2\sqrt{c}},
+\]
+!et
+and finally 
+!bt
+\[
+\frac{\partial f}{\partial x} = \frac{\partial f}{\partial a} \frac{\partial a}{\partial x}  = \frac{x(1+\exp{a})}{\sqrt{c}},
+\]
+!et
+which is just
+!bt
+\[
+\frac{\partial f}{\partial x} = \frac{x(1+b)}{d},
+\]
+!et
+and requires only three operations if we can reuse all intermediate variables.
+
+
+!split
+===== In general not this simple =====
+
+In general, see the generalization below, unless we can obtain simple
+analytical expressions which we can simplify further, the final
+implementation of automatic differentiation involves repeated
+calculations (and thereby operations) of derivatives of elementary
+functions.
+
+!split
+===== Automatic differentiation =====
+
+We can make this example more formal. Automatic differentiation is a
+formalization of the previous example (see graph).
+
+We define $\bm{x}\in x_1,\dots, x_l$ input variables to a given function $f(\bm{x})$ and $x_{l+1},\dots, x_L$ intermediate variables.
+
+In the above example we have only one input variable, $l=1$ and four intermediate variables, that is
+!bt
+\[
+\begin{bmatrix} x_1=x & x_2 = x^2=a & x_3 =\exp{a}= b & x_4=c=a+b & x_5 = \sqrt{c}=d \end{bmatrix}.
+\]
+!et
+
+Furthemore, for $i=l+1, \dots, L$ (here $i=2,3,4,5$ and $f=x_L=d$), we
+define the elementary functions $g_i(x_{Pa(x_i)})$ where $x_{Pa(x_i)}$ are the parent nodes of the variable $x_i$.
+
+In our case, we have for example for $x_3=g_3(x_{Pa(x_i)})=\exp{a}$, that $g_3=\exp{()}$ and $x_{Pa(x_3)}=a$.
+
+!split
+===== Chain rule =====
+
+We can now compute the gradients by back-propagating the derivatives using the chain rule.
+We have defined
+!bt
+\[
+\frac{\partial f}{\partial x_L} = 1,
+\]
+!et
+which allows us to find the derivatives of the various variables $x_i$ as
+!bt
+\[
+\frac{\partial f}{\partial x_i} = \sum_{x_j:x_i\in Pa(x_j)}\frac{\partial f}{\partial x_j} \frac{\partial x_j}{\partial x_i}=\sum_{x_j:x_i\in Pa(x_j)}\frac{\partial f}{\partial x_j} \frac{\partial g_j}{\partial x_i}.
+\]
+!et
+
+Whenever we have a function which can be expressed as a computation
+graph and the various functions can be expressed in terms of
+elementary functions that are differentiable, then automatic
+differentiation works.  The functions may not need to be elementary
+functions, they could also be computer programs, although not all
+programs can be automatically differentiated.
+
+!split
+===== First network example, simple percepetron with one input =====
+
+As yet another example we define now a simple perceptron model with
+all quantities given by scalars. We consider only one input variable
+$x$ and one target value $y$.  We define an activation function
+$\sigma_1$ which takes as input
+
+!bt
+\[
+z_1 = w_1x+b_1,
+\]
+!et
+where $w_1$ is the weight and $b_1$ is the bias. These are the
+parameters we want to optimize.  The output is $a_1=\sigma(z_1)$ (see
+graph from whiteboard notes). This output is then fed into the
+_cost/loss_ function, which we here for the sake of simplicity just
+define as the squared error
+
+!bt
+\[
+C(x;w_1,b_1)=\frac{1}{2}(a_1-y)^2.
+\]
+!et
+
+!split
+===== Layout of a simple neural network with no hidden layer  =====
+
+FIGURE: [figures/simplenn1.png, width=900 frac=1.0]
+
+
+
+!split
+===== Optimizing the parameters =====
+
+In setting up the feed forward and back propagation parts of the
+algorithm, we need now the derivative of the various variables we want
+to train.
+
+We need
+!bt
+\[
+\frac{\partial C}{\partial w_1} \hspace{0.1cm}\mathrm{and}\hspace{0.1cm}\frac{\partial C}{\partial b_1}. 
+\]
+!et
+
+Using the chain rule we find 
+!bt
+\[
+\frac{\partial C}{\partial w_1}=\frac{\partial C}{\partial a_1}\frac{\partial a_1}{\partial z_1}\frac{\partial z_1}{\partial w_1}=(a_1-y)\sigma_1'x,
+\]
+!et
+and
+!bt
+\[
+\frac{\partial C}{\partial b_1}=\frac{\partial C}{\partial a_1}\frac{\partial a_1}{\partial z_1}\frac{\partial z_1}{\partial b_1}=(a_1-y)\sigma_1',
+\]
+!et
+which we later will just define as
+!bt
+\[
+\frac{\partial C}{\partial a_1}\frac{\partial a_1}{\partial z_1}=\delta_1.
+\]
+!et
+
+
+!split
+===== Adding a hidden layer =====
+
+We change our simple model to (see graph)
+a network with just one hidden layer but with scalar variables only.
+
+Our output variable changes to $a_2$ and $a_1$ is now the output from the hidden node and $a_0=x$.
+We have then
+!bt
+\[
+z_1 = w_1a_0+b_1 \hspace{0.1cm} \wedge a_1 = \sigma_1(z_1),
+\]
+!et
+!bt
+\[
+z_2 = w_2a_1+b_2 \hspace{0.1cm} \wedge a_2 = \sigma_2(z_2),
+\]
+!et
+and the cost function
+!bt
+\[
+C(x;\bm{\Theta})=\frac{1}{2}(a_2-y)^2,
+\]
+!et
+with $\bm{\Theta}=[w_1,w_2,b_1,b_2]$.
+
+
+!split
+===== Layout of a simple neural network with one hidden layer  =====
+
+FIGURE: [figures/simplenn2.png, width=900 frac=1.0]
+
+
+
+!split
+===== The derivatives =====
+
+The derivatives are now, using the chain rule again
+
+!bt
+\[
+\frac{\partial C}{\partial w_2}=\frac{\partial C}{\partial a_2}\frac{\partial a_2}{\partial z_2}\frac{\partial z_2}{\partial w_2}=(a_2-y)\sigma_2'a_1=\delta_2a_1,
+\]
+!et
+!bt
+\[
+\frac{\partial C}{\partial b_2}=\frac{\partial C}{\partial a_2}\frac{\partial a_2}{\partial z_2}\frac{\partial z_2}{\partial b_2}=(a_2-y)\sigma_2'=\delta_2,
+\]
+!et
+!bt
+\[
+\frac{\partial C}{\partial w_1}=\frac{\partial C}{\partial a_2}\frac{\partial a_2}{\partial z_2}\frac{\partial z_2}{\partial a_1}\frac{\partial a_1}{\partial z_1}\frac{\partial z_1}{\partial w_1}=(a_2-y)\sigma_2'a_1\sigma_1'a_0,
+\]
+!et
+!bt
+\[
+\frac{\partial C}{\partial b_1}=\frac{\partial C}{\partial a_2}\frac{\partial a_2}{\partial z_2}\frac{\partial z_2}{\partial a_1}\frac{\partial a_1}{\partial z_1}\frac{\partial z_1}{\partial b_1}=(a_2-y)\sigma_2'\sigma_1'=\delta_1.
+\]
+!et
+
+Can you generalize this to more than one hidden layer?
+
+
+!split
+===== Important observations =====
+
+!bblock
+From the above equations we see that the derivatives of the activation
+functions play a central role. If they vanish, the training may
+stop. This is called the vanishing gradient problem, see discussions below. If they become
+large, the parameters $w_i$ and $b_i$ may simply go to infinity. This
+is referenced as  the exploding gradient problem.
+!eblock
+
+
+
+
+!split
+===== The training =====
+
+The training of the parameters is done through various gradient descent approximations with
+
+!bt
+\[
+w_{i}\leftarrow w_{i}- \eta \delta_i a_{i-1},
+\]
+!et
+and
+!bt
+\[
+b_i \leftarrow b_i-\eta \delta_i,
+\]
+!et
+with $\eta$ is the learning rate.
+
+One iteration consists of one feed forward step and one back-propagation step. Each back-propagation step does one update of the parameters $\bm{\Theta}$.
+
+For the first hidden layer $a_{i-1}=a_0=x$ for this simple model.
+
+