Convolutional Neural Nets in Chemical Engineering:
Foundations, Computations, and Applications

The convolution of scalar continuous functions $u:\mathbb{R}\mapsto\mathbb{R}$ and $v:\mathbb{R}\mapsto\mathbb{R}$ is denoted as $\psi=u*v$ and their cross-correlation is denoted as $\phi=u\star v$. These operations are given by:

We will refer to function $u$ as the convolution operator and to $v$ as the input signal. The output of the convolution operation is a scalar continuous function $\psi:\mathbb{R}\mapsto\mathbb{R}$ that we refer to as the convolved signal. The output of the cross-correlation operation is also a continuous function $\phi:\mathbb{R}\mapsto\mathbb{R}$ that we refer to as the cross-correlated signal.

The convolution and cross-correlation are applied to the signal by spanning the domain $x\in(-\infty,\infty)$; we can see that convolution is applied by looking backward, while cross-correlation is applied by looking forward. One can show that cross-correlation is equivalent to convolution that uses a rotated (flipped) operator $u$ by 180 degrees; in other words, if we define the rotated operator as $-u$ then $u*v=(-u)\star v$. We thus have that convolution and cross-correlation are mirror operations; consequently, convolution and cross-correlation are often both referred to as convolution. In what follows we use the term convolution and cross-correlation interchangeably; we highlight one operation over the other when appropriate. Modern CNN packages such as PyTorch [29] use cross-correlation in their implementation (this is more intuitive and easier to implement).

One can think of convolution and cross-correlation as weighting operations (with weighting function defined by the operator). The weighting function $u$ can be designed in different ways to perform different operations to the signal $v$ such as averaging, smoothing, differentiation, and pattern recognition. In Figure 1 we illustrate the application of a Gaussian operator $u$ to a noisy rectangular signal $v$. We note that the convolved and cross-correlated signals are identical (because the operator is symmetric and thus $u=-u$). We also note that the output signals are smooth versions of the input signal; this is because a Gaussian operator has the effect of extracting frequencies from the signal (a fact that can be established from Fourier analysis). Specifically, we recall that the Fourier transform of the convolved signal $\psi$ satisfies:

where $\mathcal{F}\{w\}$ is the Fourier transform of scalar function $w$. Here, the product $\mathcal{F}\{u\}\cdot\mathcal{F}\{v\}$ has the effect of preserving the frequencies in the signal $v$ that are also present in the operator $u$. In other words, convolution acts as a filter of frequencies; as such, one often refers to operators as filters. By comparing the output signals of convolution and cross-correlation obtained with the triangular operator, we can confirm that one is the mirror version of the other. From Figure 1, we also see that the output signals are smooth variants of the input signal; this is because the triangular operator also extracts frequencies from the input signal (but the extraction is not as clean as that obtained with the Gaussian operator). This highlights the fact that different operators have different frequency content (known as the frequency spectrum).

Convolution and cross-correlation are implemented computationally by using discrete representations; these signals are represented by column vectors $u\in\mathbb{R}^{n}$ and $v\in\mathbb{R}^{n}$ with entries denoted as $u[x]$, $v[x]$, $x=-N,...,N$ (note that $n=2\cdot N+1$). Here, we note that the input and operator are defined over a 1D grid. The discrete convolution results in vector $\psi\in\mathbb{R}^{n}$ with entries given by:

Here, we use $\{-N,N\}$ to denote the integer sequence $-N,-N+1,...,-1,0,1,...,N-1,N$. We note that computing $\psi\left[x\right]$ at a given location $x$ requires $n$ operations and thus computing the entire output signal (spanning $x\in\{-N,N\}$) requires $n^{2}$ operations; consequently, convolution becomes expensive when the signal and operator are long vectors; as such, the operator $u$ is often a vector of low dimension (compared to the dimension of the signal $v$). Specifically, consider the operator $u\in\mathbb{R}^{n_{u}}$ with entries $u[x],\,x=-N_{u},...,N_{u}$ (thus $n_{u}=2\cdot N_{u}+1$) and assume $n_{u}\ll n$ (typically $n_{u}=3$). The convolution is given by:

This reduces the number of operations needed from $n^{2}$ to $n\cdot n_{u}$.

The convolved signal is obtained by using a moving window starting at the boundary $x=-N$ and by moving forward as $x+1$ until reaching $x=N$. We note that the convolution is not properly defined close to the boundaries (because the window lies outside the domain of the signal). This situation can be remedied by starting the convolution at an inner entry of $v$ such that the full window is within the signal (this gives a proper convolution). However, this approach has the effect of returning a signal $\psi$ that is smaller than the original signal $v$. This issue can be overcome by artificially padding the signal by adding zero entries (i.e., adding ghost entries to the signal at the boundaries). This is an improper convolution but returns a signal $\psi$ that has the same dimension as $v$ (this is often convenient for analysis and computation). Figure 2 illustrates the difference between convolutions with and without padding. We also highlight that it is possible to advance the moving window by multiple entries as $x+s$ (here, $s\in\mathbb{Z}_{+}$ is known as the stride). This has the effect of reducing the number of operations (by skipping some entries in the signal) but returns a signal that is smaller than the original one.

By defining signals $u\in\mathbb{R}^{n_{u}}$ and $v\in\mathbb{R}^{n_{v}}$ with entries $u[x],\ x=1,...,n_{u}$ and $v[x],\ x=1,...,n_{v}$, one can also express the valid convolution operation in an asymmetric form as:

In CNNs, one often applies multiple operators to a single input signal or one analyzes input signals that have multiple channels. This gives rise to the concepts of single-input single-output (SISO), single-input multi-output (SIMO), and multi-input multi-output (MIMO) convolutions.

In SISO convolution, we take a 1-channel input vector and output a vector (a 1-channel signal), as in (5). In SIMO convolution, one uses a single input $v\in\mathbb{R}^{n_{v}}$ and a collection of $q$ operators $\mathcal{U}_{\left(j\right)}\in\mathbb{R}^{n_{u}}$ with $j\in\{1,q\}$. Here, every element $\mathcal{U}_{\left(j\right)}$ is an $n_{u}$-dimensional vector and we represent the collection as the object $\mathcal{U}$. SIMO convolution yields a collection of convolved signals $\Psi_{\left(j\right)}\in\mathbb{R}^{n_{v}-n_{u}+1}$ with $j\in\{1,q\}$ and with entries given by:

In MIMO convolution, we consider a multi-channel input given by the collection $\mathcal{V}_{\left(i\right)}\in\mathbb{R}^{n_{v}},\,i\in\{1,p\}$ ($p$ is the number of channels); the collection is represented as the object $\mathcal{V}$. We also consider a collection of operators $\mathcal{U}_{(i,j)}\in\mathbb{R}^{n_{u}}$ with $i\in\{1,p\}$ and $j\in\{1,q\}$; in other words, we have $q$ operators per input channel $i\in\{1,p\}$ and we represent the collection as the object $\mathcal{U}$. MIMO convolution yields the collection of convolved signals $\Psi_{\left(j\right)}\in\mathbb{R}^{n_{v}-n_{u}+1}$ with $j\in\{1,q\}$ and with entries given by:

We see that, in MIMO convolution, we add the contribution of all channels $i\in\{1,p\}$ to obtain the output object $\Psi$, which contains $j\in\{1,q\}$ channels given by vectors $\Psi_{(i)}$ of dimension $(n_{v}-n_{u}+1)$. Channel combination loses information but saves computer memory.

Convolution and cross-correlation operations in 2D are analogous to those in 1D; the convolution and cross-correlation of a continuous operator $u:\mathbb{R}^{2}\mapsto\mathbb{R}$ and an input signal $v:\mathbb{R}^{2}\mapsto\mathbb{R}$ are denoted as $\psi=u*v$ and $\phi=u\star v$ and are given by:

As in the 1D case, the terms convolution and cross-correlation are used interchangeably; here, we will use the cross-correlation form (typically used in CNNs).

In the discrete case, the input signal and the convolution operator are matrices $U\in\mathbb{R}^{n_{U}\times n_{U}}$ and $V\in\mathbb{R}^{n_{V}\times n_{V}}$ with entries $U\left[x_{1},x_{2}\right],\;x_{1},x_{2}\in\{-N_{U},N_{U}\}$ and $V\left[x_{1},x_{2}\right],\;x_{1},x_{2}\in\{-N_{V},N_{V}\}$ (thus $n_{U}=2\cdot N_{U}+1$ and $n_{V}=2\cdot N_{V}+1$). For simplicity, we assume that these are square matrices and we note that the input and operator are defined over a 2D grid. The convolution of the input and operator matrices results in a matrix $\Psi=U*V$ with entries:

In Figure 3 we illustrate 2D convolutions with and without padding. The 2D convolution (in valid and asymmetric form) can be expressed as:

where $x_{1}\in\{1,n_{V}-n_{U}+1\}$, $x_{2}\in\{1,n_{V}-n_{U}+1\}$ and $\Psi\in\mathbb{R}^{\left(n_{V}-n_{U}+1\right)\times\left(n_{V}-n_{U}+1\right)}$.

The SISO convolution of an input $V\in\mathbb{R}^{n_{V}\times n_{V}}$ and an operator $U\in\mathbb{R}^{n_{U}\times n_{U}}$ is given by (10) and outputs a matrix $\Psi\in\mathbb{R}^{\left(n_{V}-n_{U}+1\right)\times\left(n_{V}-n_{U}+1\right)}$. In SIMO convolution, we are given a collection of operators $\mathcal{U}_{\left(j\right)}\in\mathbb{R}^{n_{\mathcal{U}}\times n_{\mathcal{U}}}$ with $j\in\{1,q\}$. A convolved matrix $\Psi_{\left(j\right)}\in\mathbb{R}^{\left(n_{V}-n_{\mathcal{U}}+1\right)\times\left(n_{V}-n_{\mathcal{U}}+1\right)}$ is obtained by applying the $j$-th operator $\mathcal{U}_{\left(j\right)}$ to the input $V$:

for $j\in\{1,q\},\;x_{1},x_{2}\in\{1,n_{V}-n_{\mathcal{U}}+1\}$. The collection of convolved matrices $\Psi_{\left(j\right)}\in\mathbb{R}^{\left(n_{V}-n_{\mathcal{U}}+1\right)\times\left(n_{V}-n_{\mathcal{U}}+1\right)},\;\ j\in\{1,q\}$ is represented as object $\Psi$.

In MIMO convolution, we are given a $p$-channel input collection $\mathcal{V}_{\left(i\right)}\in\mathbb{R}^{n_{\mathcal{V}}\times n_{\mathcal{V}}}$ with $i\in\{1,p\}$ and represented as $\mathcal{V}$. We convolve this input with the operator object $\mathcal{U}$, which is a collection $\mathcal{U}_{(i,j)}\in\mathbb{R}^{n_{\mathcal{U}}\times n_{\mathcal{U}}},\;i\in\{1,p\},\;j\in\{1,q\}$. This results in an object $\Psi$ given by the collection $\Psi_{(j)}\in\mathbb{R}^{\left(n_{\mathcal{V}}-n_{\mathcal{U}}+1\right)\times\left(n_{\mathcal{V}}-n_{\mathcal{U}}+1\right)}$ for $j\in\{1,q\}$ and with entries given by:

for $j\in\{1,q\},\;x_{1}\in\{1,n_{\mathcal{V}}-n_{\mathcal{U}}+1\}$, and $x_{2}\in\{1,n_{\mathcal{V}}-n_{\mathcal{U}}+1\}$. For convenience, the input object is often represented as a 3D tensor $\mathcal{V}\in\mathbb{R}^{n_{\mathcal{V}}\times n_{\mathcal{V}}\times p}$, the convolution operator is represented as the 4D tensor $\mathcal{U}\in\mathbb{R}^{n_{\mathcal{U}}\times n_{\mathcal{U}}\times p\times q}$, and the output signal is represented as the 3D tensor $\Psi\in\mathbb{R}^{\left(n_{\mathcal{V}}-n_{\mathcal{U}}+1\right)\times\left(n_{\mathcal{V}}-n_{\mathcal{U}}+1\right)\times q}$. We note that, if $p=1$ and $q=1$ (1-channel inputs and channels), these tensors become matrices. Tensors are high-dimensional quantities that require significant computer memory to store and significant power to process.

MIMO convolutions in 2D are typically used to process RGB images (3-channel input), as shown in Figure 4. Here, the RGB image is the object $\mathcal{V}$ and each input channel $\mathcal{V}_{\left(i\right)}$ is a matrix; each of these matrices is convolved with an operator $\mathcal{U}_{\left(i,1\right)}$. Here we assume $q=1$ (one operator per channel). The collection of convolved matrices $\Psi_{\left(i,1\right)}$ are combined to obtain a single matrix $\Psi_{\left(1\right)}$. If we consider a collection of operators $\mathcal{U}_{\left(i,j\right)}$ with $i\in\{1,p\}$, $j\in\{1,q\}$ and $q>1$, the output of MIMO convolution returns the collection of matrices $\Psi_{\left(j\right)}$ with $j\in\{1,q\}$, which is assembled in the tensor $\Psi$. Convolution with multiple operators allows for the extraction of different types of features from different channels.

We highlight that the use of channels is not restricted to images; specifically, channels can be used to input data of different variables in a grid (e.g., temperature, concentration, density). As such, channels provide a flexible framework to express multivariate data.

We also highlight that a grayscale image is a 1-channel input matrix (the RGB channels are combined in a single channel). In a grayscale image, every pixel (an entry in the associated matrix) has a certain light intensity value; whiter pixels have higher intensity and darker pixels have a lower intensity (or the other way around). The resolution of the image is dictated by the number of pixels and thus dictates the size of the matrix; the size of the matrix dictates the amount of memory needed for storage and computations needed for convolution. It is also important to emphasize that any matrix can be visualized as a grayscale image (and any grayscale image has an underlying matrix). This duality is important because visualizing large matrices (e.g., to identify patterns) is difficult if one directly inspects the numerical values; as such, a convenient approach to analyze patterns in large matrices consists of visualizing them as images.

Finally, we highlight that convolution operators can be applied to any grid data object in higher dimensions (e.g., 3D and 4D) in a similar manner. Convolutions in 3D can be used to process video data (each time frame is an image). However, 3D data objects can also be used to represent data distributed over 3D Euclidean spaces (e.g., density or flow in a 3D domain). However, the complexity of convolution operations in 3D (and higher dimensions) is substantial. In the discussion that follows we focus our attention to 2D convolutions; in Section 6 we illustrate the use of 3D convolutions in a practical application.

An input sample of the CNN is the tensor $\mathcal{V}\in\mathbb{R}^{n_{\mathcal{V}}\times n_{\mathcal{V}}\times p}$. The convolution block uses the operator $\mathcal{U}\in\mathbb{R}^{n_{\mathcal{U}}\times n_{\mathcal{U}}\times p\times q}$ to conduct a MIMO convolution. The output of this convolution is the tensor $\Psi\in\mathbb{R}^{n_{\Psi}\times n_{\Psi}\times q}$ with entries given by:

where $x_{1}\in\{1,n_{\Psi}\}$, $x_{2}\in\{1,n_{\Psi}\}$, and $j\in\{1,q\}$. A bias parameter $b_{c}[j]\in\mathbb{R}$ is added after the convolution operation; the bias helps adjust the magnitude of the convolved signal. To enable compact notation, we define the convolution block using the mapping $\Psi=f_{c}(\mathcal{V};\mathcal{U},b_{c})$; here, we note that the mapping depends on the parameters $\mathcal{U}$ (the operators) and $b_{c}$ (the bias). An example of a convolution block with a 3-channel input and 2-channel operator is shown in Figure 7.

The convolution block outputs the signal $\Psi$; this is passed through an activation block given by the mapping $\mathcal{A}=h_{c}(\Psi)$ with $\mathcal{A}\in\mathbb{R}^{n_{\Psi}\times n_{\Psi}\times q}$. Here, we define an activation mapping of the form $h_{c}:\mathbb{R}^{n_{\Psi}\times n_{\Psi}\times q}\mapsto\mathbb{R}^{n_{\Psi}\times n_{\Psi}\times q}$. The activation of $\Psi$ is conducted element-wise as:

where $\alpha:\mathbb{R}\mapsto\mathbb{R}$ is an activation function (a scalar function). Typical activation functions include the sigmoid, hyperbolic tangent (tanh), and Rectified Linear Unit (ReLU):

Figures 8 and 9 illustrate the transformation induced by the activation functions. These functions act as basis functions that, when combined in the CNN, enable capturing nonlinear behavior. A common problem with sigmoid and tanh functions is that they exhibit the so-called vanishing gradient effect [11]. Specifically, when the input values are large or small in magnitude, the gradient of the sigmoid and tanh functions is small (flat at both ends) and this makes the activation output insensitive to changes in the input. Furthermore, both the sigmoid and tanh functions are sensitive to the change in input when the output is close to 1/2 and zero, respectively. The ReLU function is commonly used to avoid vanishing gradient effects and to increase sensitivity [26]. This activation function outputs a value of zero when the input is less than or equal to zero, but outputs the input value itself when the input is greater than zero. The function is linear when the input is greater than zero, which makes the CNN easier to optimize with gradient-based methods [11]. However, ReLU is not continuously differentiable when the input is zero; in practice, CNN implementations assume that the gradient is zero when the input is zero.

The pooling block is a transformation that reduces the size of the output obtained by convolution and subsequent activation. This block also seeks to make the representation approximately invariant to small translation of the inputs [11]. Max-pooling and average-pooling are the most common pooling operations (here we focus on max-pooling). The pooling operation $\mathcal{P}=f_{p}\left(\mathcal{A}\right)$ can be expressed as a mapping $f_{p}:\mathbb{R}^{n_{\Psi}\times n_{\Psi}\times q}\mapsto\mathbb{R}^{\left\lfloor n_{\Psi}/n_{p}\right\rfloor\times\left\lfloor n_{\Psi}/n_{p}\right\rfloor\times q}$ and delivers a tensor $\mathcal{P}\in\mathbb{R}^{\left\lfloor n_{\Psi}/n_{p}\right\rfloor\times\left\lfloor n_{\Psi}/n_{p}\right\rfloor\times q}$. To simplify the discussion, we denote $n_{\mathcal{P}}=\left\lfloor n_{\Psi}/n_{p}\right\rfloor$; the max-pooling operation with $i\in\{1,q\}$ is defined as:

where $x_{1}\in\left\{1,n_{\mathcal{P}}\right\}$, $x_{2}\in\left\{1,n_{\mathcal{P}}\right\}$. In Figure 10, we illustrate the max-pooling and averaging pooling operations on a 1-channel input.

Operations of convolution, activation, and pooling constitute a convolution unit and deliver an output object $\mathcal{P}=f_{p}(h_{c}(f_{c}(\mathcal{V}))$. This object can be fed into another convolution unit to obtain a new output object $\mathcal{P}=f_{p}(h_{c}(f_{c}(\mathcal{P}))$ and this recursion can be performed over multiple units. The recursion has the effect of highlighting different features of the input data object $\mathcal{V}$ (e.g., that capture local and global patterns). In our discussion we consider a single convolution unit.

The convolution, activation, and pooling blocks deliver a tensor $\mathcal{P}\in\mathbb{R}^{n_{\mathcal{P}}\times n_{\mathcal{P}}\times p}$ that is flattened to a vector $v\in\mathbb{R}^{n_{\mathcal{P}}\cdot n_{\mathcal{P}}\cdot p}$ (with some abuse of notation where $v$ was originally defined for input signal in 1D); this vector is fed into a fully-connected (dense) block that performs the final prediction. The vector $v$ is typically known as the feature vector. The flattening block is represented as the mapping $f_{f}:\mathbb{R}^{n_{\mathcal{P}}\times n_{\mathcal{P}}\times p}\mapsto\mathbb{R}^{n_{v}}$ with $n_{v}=n_{\mathcal{P}}\cdot n_{\mathcal{P}}\cdot p$ that outputs $v=f_{f}\left(\mathcal{P}\right)$. Note that this block is simply a vectorization of a tensor.

The feature vector $v\in\mathbb{R}^{n_{v}}$ is input into a prediction block that delivers the final prediction $\hat{y}\in\mathbb{R}$. This block is typically a fully-connected (dense) neural net with multiple weighting and activation units (here we only consider a single unit). The prediction block first mixes the elements of the feature vector as:

where $w\in\mathbb{R}^{n_{v}}$ is the weighting (parameter) vector and $b_{d}\in\mathbb{R}$ is the bias parameter. The scalar $d\in\mathbb{R}$ is normally known as the evidence. We denote the weighting operation using the mapping $f_{d}:\mathbb{R}^{n_{v}}\mapsto\mathbb{R}$ and thus $d=f_{d}(v;w,b_{d})$. In Figure 11, we illustrate this weighting operation. An activation function is applied to the evidence $d$ to obtain the final prediction $\hat{y}=h_{d}(d)$ with $h_{d}:\mathbb{R}\mapsto\mathbb{R}$. This involves an activation (e.g., with a sigmoidal or ReLU function). We thus have that the prediction block has the form $\hat{y}=h_{d}(f_{d}(v))$. One can build a multi-layer fully-connected network by using a recursion of weighting and activation steps (here we consider one layer to simplify the exposition). Here, we assume that the CNN delivers a single (scalar) output $\hat{y}$. In practice, however, a CNN can also predict multiple outputs; in such a case, the weight parameter becomes a matrix and the bias is a vector. We also highlight that the predicted output can be an integer value (to perform classification tasks) or a continuous value (to perform regression tasks).

To explain the process of forward propagation, we consider the CNN architecture shown in Figure 6 with an input $\mathcal{V}\in\mathbb{R}^{n_{\mathcal{V}}\times n_{\mathcal{V}}\times p}$ and an output (prediction) $\hat{y}\in\mathbb{R}$. All parameters $\mathcal{U}\in\mathbb{R}^{n_{\mathcal{U}}\times n_{\mathcal{U}}\times p\times q}$, $b_{c}\in\mathbb{R}^{q}$, $w\in\mathbb{R}^{n_{v}}$ and $b_{d}\in\mathbb{R}$ are flattened and concatenated to form a parameter vector $\theta\in\mathbb{R}^{n_{\theta}}$, where $n_{\theta}=n_{\mathcal{U}}\cdot n_{\mathcal{U}}\cdot p\cdot q+q+n_{\mathcal{V}}+1$. We can express the entire sequence of operations carried out in the CNN as a composite mapping $F:\mathbb{R}^{n_{\mathcal{V}}\times n_{\mathcal{V}}\times p}\mapsto\mathbb{R}$ and thus $\hat{y}=F\left(\mathcal{V};\theta\right)$. We call this mapping the forward propagation mapping (shown in Figure 6) and note that this can be decomposed sequentially (via lifting) as:

This reveals the nested nature and the order of the CNN transformations on the input data. Specifically, we see that one conducts convolution, activation, pooling, flattening, and prediction.

To perform training, we collect a set of input samples $\mathcal{V}^{(i)}\in\mathbb{R}^{n_{\mathcal{V}}\times n_{\mathcal{V}}\times p}$ and corresponding output labels $y[i]\in\mathbb{R}$ with sample index $i\in\{1,n\}$. Our aim is to compare the true output labels $y[i]$ to the CNN predictions $\hat{y}[i]=F\left(\mathcal{V}^{(i)};\theta\right)$. To do so, we define the loss function $\mathcal{L}:\mathbb{R}^{n}\to\mathbb{R}$ that we aim to minimize. For classification tasks, the output $\hat{y}[i]$ is binary and a common loss function is the cross-entropy [17]:

For regression tasks, we have that the output $\hat{y}^{(i)}$ is continuous and a common loss function is the mean squared error:

In compact form, we can write the loss function as $\mathcal{L}\left(F\left(\mathcal{V};\theta\right)\right)=\mathcal{L}\left(\theta\right)$. Here, $\mathcal{V}$ contains all the input samples $\mathcal{V}^{(i)}$ for $i\in\{1,n\}$. The loss function can be decomposed into individual samples and thus we can express it as:

where $\mathcal{L}^{(i)}(\theta)$ is the contribution associated with the $i$-th data sample.

The parameters of the CNN are obtained by finding a solution to the optimization problem:

The loss function $\mathcal{L}$ embeds the forward propagation mapping $F$ (which is a highly nonlinear mapping); as such, the loss function might contain many local minima (and such minima tend to be flat [12]). We recall that appearance of flat minima is often the manifestation of model overparameterization; this is due to the fact that the parameter vector $\theta$ can contain hundreds of thousands to millions of values. As such, training of CNNs need to be equipped with appropriate validation procedures (to ensure that the model generalizes well). Most optimization algorithms used for CNNs are gradient-based (first-order); this is because second-order optimization methods (e.g., Newton-based) involve calculating the Hessian of the loss function $\mathcal{L}$ with respect to parameter vector $\theta$, which is computationally expensive (or intractable). Quasi-Newton methods (e.g., limited-memory BFGS) are also often used to train CNNs.

Gradient descent (GD) is the most commonly used algorithm to train CNNs. This algorithm updates the parameters as:

where $\eta\in\mathbb{R}_{+}$ is the learning rate (step length) and $\Delta\theta\in\mathbb{R}^{n_{\theta}}$ is the step direction. In GD, the step direction is set to $\Delta\theta=\nabla_{\theta}\mathcal{L}(\theta)$ (the gradient of the loss). Since the loss function is a sum over samples, the gradient can be decomposed as:

Although GD is easy to compute and implement, $\theta$ is not updated until the gradient of the entire dataset is calculated (which can contain millions of data samples). In other words, the parameters are not updated until all the gradient components $\nabla_{\theta}\mathcal{L}^{(i)}(\theta)$ are available. This leads to load-imbalancing issues and limits the parallel scalability of the algorithm.

Stochastic gradient descent (SGD) is a variant of GD that updates $\theta$ more frequently (with gradients from a partial number of samples). Specifically, $\theta$ changes after calculating the loss for each sample. SGD updates the parameters by using the step $\Delta\theta=\nabla_{\theta}\mathcal{L}^{(i)}(\theta)$, where sample $i$ is selected at random. Since $\theta$ is updated frequently, SGD requires less memory and has better parallel scalability [20]. Moreover, it has been shown empirically that this algorithm has the tendency to escape local minima (it explores the parameter space better). However, if the training sample has high variance, SGD will converge slowly.

Mini-batch GD is an enhancement of SGD; here, the entire dataset is divided into $b$ batches and $\theta$ is updated after calculating the gradient of each batch. This results in a step direction of the form:

where $\mathcal{B}^{(i)}$ is a set of sample corresponding to batch $i\in\{1,b\}$ and the entries of the batches are selected at random. Mini-batch GD updates the model parameters frequently and has faster convergence in the presence of high variance [41].

Backward propagation (backpropagation) seeks to compute elements of the gradient of the loss $\nabla_{\theta}\mathcal{L}^{(i)}(\theta)$ by recursive use of the chain rule. This approach exploits the nested nature of the forward propagation mapping:

This mapping can be expressed in backward form as:

An illustration of this backward sequence is presented in Figure 6. Exploiting this structure is essential to enable scalability (particularly for saving memory). To facilitate the explanation, we consider the squared error loss defined for a single sample (generalization to multiple samples is trivial due to the additive nature of the loss):

Our goal is to compute the gradient (the parameter step direction) $\Delta\theta=\nabla\mathcal{L}(\theta)$; here, we have that parameter $\theta=(w,b_{d},\mathcal{U},b_{c})$.