In our previous post, we discussed how to compute the output of a network using a given set of parameters (weights and biases) along with input data. This process is known as the forward pass or forward propagation.

To determine the set of parameters that best models the data, we must solve an optimization problem. This post will explain how to approach this problem and apply the method to a spiral data classification case.

The learning problem

Our previous post observed the network’s capability once it was trained. Let’s revisit the multi-layer perceptron (MLP) model, which has an input/output relationship utilizing, for example, two hidden layers:

For a given input x, the output provides the model prediction, which could correspond to one of the three colors in the spiral dataset. To ensure the model predicts well, we need to evaluate its accuracy using a cost function. The cost function measures the error between the model’s prediction and the actual training data. Let us say we have N training data.

A simple choice for the cost function is the sum of the squared errors for each training data, labeled by 1 ≤ i ≤ N

This cost function calculates the least-squares distance between the predicted and true values, emphasizing the dependency on all the model’s parameters — weights and biases — that need to be optimized.

A high value of the cost function indicates poor model predictions, while a low value suggests that the model is likely predicting well. The goal is to find the set of weights and biases that minimize this cost.

Note that for the spiral dataset, we will use a different cost function, which we will discuss later in the post.

A recipe to minimize the cost function

Minimizing the cost function is a numerical optimization problem and is a complex topic. We describe a simple way to do so, which is standard in most deep learning problems.

The idea is to follow an iterative process: we randomly choose a set of starting parameters (weights and biases) of the model, and update them at each iteration to gradually reduce the cost function. The key question is: how can we update these parameters to decrease the cost?

To answer this, we need to understand how sensitive the cost function is to changes in the parameters. The gradient of the cost function captures this sensitivity. The gradient is a vector that provides the direction and rate of the steepest increase in the cost function. Each component of this vector indicates how much the cost would change if a particular parameter were adjusted.

For a model with parameters (w, b), the gradient provides the information needed to update these parameters in a way that reduces the cost.

We now have a recipe to minimize the cost function as an iterative process:

1. At each iteration, start from an initial guess for the parameters.
2. Compute the gradient of the cost function with respect to each parameter.
3. Update the parameters in the direction that decreases the cost.

This method is known as gradient descent. At each iteration, the idea is to update the next bias and weights as follows

where η is a positive number known as the learning rate or step size, which controls the size of each update. While we fix this quantity for simplicity, more sophisticated methods adjust the learning rate dynamically.

Analogy: Imagine you’re in the dark on a steep mountain and want to descend to a village below. You would use a stick to probe the ground and step in the direction where the slope decreases the most. You can interpret the picture above as being with a bunch of friends, starting from different places, looking for a way down to the village. Similarly, in gradient descent, you update parameters by following the direction of the steepest descent to minimize the cost function.

In deep learning, the cost function landscape can be extremely complex due to the model’s complexity and the large number of parameters. As you can see in the illustration there is no guarantee of finding the global minimum, gradient descent is quite effective in practice and forms the basis for more advanced optimization techniques. Because computing the gradient can be computationally expensive, deep learning methods often use stochastic gradient descent (SGD), where the gradient is estimated using a subset of the training data.

Computing the gradient — Automatic differentiation

The gradient of the cost function is crucial for solving deep learning problems. The next challenge consists in how to compute the gradient.

For simple models, like the spiral classification case, we can manually compute the gradient. However, for more complex models, this manual approach is impractical. Fortunately, automatic differentiation (AD) provides a powerful solution. In this post, we’ll give a brief overview of how automatic differentiation works.

Automatic differentiation follows a divide-and-conquer strategy by breaking the mathematical expression of the cost function into elementary operations (+, -, *, /, sqrt, log, exp, etc.) through a computational graph. The two key components of automatic differentiation are:

1. Differentiation Rules for Elementary Operations: These rules define how to compute derivatives for basic operations.
2. Application of the Chain Rule: The chain rule allows us to compute the derivative of a composite function by combining the derivatives of its constituent tasks.

Once we have this in mind, we can break any mathematical expression (such as given by a neural network) into a computational graph, and compute the derivatives between connected nodes for each operation.

Here’s how automatic differentiation works:

1. Construct a Computational Graph (Divide): Decompose the mathematical expression of the cost function into a graph where nodes represent operations and edges represent dependencies. For example, for the expression e = c*d=(a + b)*(b + 1), you create a graph showing how a and b contribute to e through the intermediate quantities c and d.
2. Apply Differentiation Rules: Use predefined differentiation rules for each operation between connected nodes in the graph. In our example, we need the differentiation rule for multiplication and addition. The first connection gives de/dc = d(c*d)/dc = d. The next one gives dc/da=d(a+b)/da = 1, etc. Do so and fill each edge of the graph.
3. Gradient Accumulation (Conquer): At the end of the day, we want to know how changes in a and b affect e, that is de/da and de/db. Thanks to the chain rule, we can accumulate the results from step 2 to obtain the global changes. Applied to a the chain rule gives de/da = de/dc*dc/da, which from step 2 gives d*1=b+1. Applied to b the chain rule sums from both edges of the graph and we get de/db = de/dc*dc/db + de/dd*dd/db = d*1 + c*1 = d+c = 2b + a + 1.

What we have described is known as reverse-gradient accumulation, because we start from the output e and go back to the inputs a and b. Put shortly, the process accumulates gradients through each edge by multiplying derivatives along the path, and finally add all possible path together.

This is exactly the backpropagation algorithm, it computes the gradient of the cost through reverse-gradient accumulation. The strength of the approach is that the gradient computation is exact, there is no approximation.

Automatic differentiation is implemented in most deep learning software such as Pytorch, TensorFlow, JAX, etc. With one of these tools, you can define any custom forward model and directly hit the “train” or “fit” function: this will launch an optimization algorithm such as the gradient descent, including the automatic computation of the gradient.

Training the spiral classification problem

Now that we have discussed some core concepts of deep learning, let us apply them to classify the spiral dataset.

We’ll first focus on the single-layer perceptron (SLP) and then extend the approach to a multi-layer perceptron (MLP) with one hidden layer.

Before that, we need to use a suitable cost function. For each training data, the output of the SLP or MLP is y_i, a [3x1] vector with 3 arbitrary numbers. However we would like the output to represent one of the three colors classes {red, purple, green} to compare with the exact data. This conversion is done via a softmax function.

To compute the cost, we take the natural log of the probability vector, and select the component associated to the color of the training data. Finally we average over all the training points. The mathematical expression is cumbersome, but translates as follows

This is known as the cross-entropy cost. In short, the choice of the cost function is not obvious and depends on the problem.

1. Single-Layer Perceptron (SLP)

To compute the gradient of this cost, we follow the logic of automatic differentiation. We do it by hand for the model with no hidden layer to understand how things work. We decompose the gradient of the cost according to the chain rule, compute the elementary gradients thanks to differentiation rules and multiply the results together. The output of the SLP network is y_i = W*x + b, and computations give, for the weights and biases,

Now that we have the gradient of the cost , we can implement the gradient descent iterations inside a loop to update the weights and biases. Inside the loop, we i) compute the output of the network and the cost function, ii) the derivatives of the cost, and iii) finally update the weights and biases. Here is a draft of the code

# Train a Linear Classifier

# initialize parameters randomly, 
# with here D=2 (dimension) and K=3 (number of classes)
W = 0.01 * np.random.randn(D,K) # matrix of weights
b = np.zeros((1,K)) # vector of biases

num_examples = X.shape[0] # X are our ~200 data points
for n in range(300): # gradient descent loop, 300 iterations
  # 1 - output of the single-layer perceptron given the data X
  y_i = np.dot(X, W) + b 
  # compute the class probabilities
  probs = np.exp(y_i) / np.sum(np.exp(y_i), axis=1, keepdims=True)
  # 1 - compute the average cross-entropy cost
  C_ik = -np.log(probs[range(num_examples),k]) # k gives the color of the data
  cost = np.sum(C_ik)/num_examples
  
  if n % 10 == 0:
    print ("iteration %d: loss %f" % (n, cost))
  
  # 2 - gradient computation 
  # dC/dyi
  dscores = probs
  dscores[range(num_examples),k] -= 1 # p_i - 1
  dscores /= num_examples # average over training points
  # backpropagate the gradient to parameters (W,b) + sum over training data
  dW = np.dot(X.T, dscores) # dC/dw = (p_i-1)*x^T
  db = np.sum(dscores, axis=0, keepdims=True) # dC/db = p_i - 1
  
  # 4 - gradient descent update, with learning rate eta=1
  W -= dW 
  b -= db

2. Multi-Layer Perceptron (MLP)

For the MLP model with one hidden layer, we must compute the intermediate partial derivatives for the weights and biases corresponding to the first and hidden layers (W⁰, W¹, b⁰, b¹), and also account for the derivative of the activation function in the hidden layer. Even though we can take advantage of the chain rule and reuse some computations, the derivation becomes quite cumbersome. We can now see why automatic differentiation is so useful !

We have added the contribution of the hidden layer in the forward model and the gradient computation in the code that is available here. The principle is exactly the same as for the SLP.

We let the gradient descent run to obtain our final set of parameters (w, b). We obtain a training accuracy of 50% for the single-layer model SLP and 98% for the multi-layer model MLP, and also get the boundary classification from the first part of the article :-)

That’s it! I hope you understand the basics of training a neural network. it relies on a forward model, a gradient descent algorithm and an accurate computation of the gradients of the cost function. This paves the way to study more advanced forward propagation models, such as convolutional or recurrent neural networks.

Code:

GitHub - linh-ktran/demystifying-deep-learning

Resources:

• CS231n Convolutional Neural Networks for Visual Recognition
• Calculus on Computational Graphs: Backpropagation

Deep Learning has become of utmost importance in the past years, impacting almost every engineering and scientific discipline. In this post, we aim to clarify what is hidden behind deep learning and review its core underlying principles. This post only assumes knowledge of basic calculus.

What do we mean by Deep Learning?

If we look at Wikipedia, Deep Learning is defined as a subset of machine learning methods based on neural networks with representation learning.

It is fair to think of deep learning as a specific category of machine learning models. Deep learning uses representation learning, which means these models are designed to learn features from data. A feature is simply a measurable property (a numeric value) that characterizes a dataset: it can be an age, a color, an animal type, a movie gender, etc.

So what do we mean by learning? Learning refers to finding an accurate representation, or an accurate model of the data. For example, a famous learning task consists of recognizing handwritten digits from a series of images.

In this task, we want to find a mathematical model that can distinguish between 10 digits. This is usually achieved in two steps :

1. Training: We fit a model to a given dataset. The model is a mathematical function described by a number of unknown parameters.
2. Evaluation: We test our model on unknown data.

If the model can evaluate the correct output, we can say that the machine has learned a representation of handwritten digits.

Deep learning is a subset of machine learning and uses a specific mathematical model for this task called artificial neural networks. The rest of this post is devoted to how to train/fit a basic artificial neural network to a given dataset, such as handwritten digits.

Artificial neural networks

Artificial neural networks are a mathematical model inspired by the neurons of the human brain. They relate input data x (e.g. the pixel grayscale value of a handwritten digit image) to output data y (e.g. an integer from 0 to 9) through a function f such that: y = f(x).

An artificial neural network, or neural net, refers to a specific structure of this function. We explain the concept of a multi-layer network in two building blocks.

1. The Artificial Neuron Model (Perceptron): An output value y is related to an input vector x = (x_1, x_2, …, x_n) by a weighted sum, parametrized by a vector of weights w = (w_1, w_2, …, w_n) and a bias value b. The weighted sum is then passed through an activation function g. The activation function, which adds nonlinear behavior, is a choice of the model. For this post, we do not explain how to choose it.

The perceptron encodes a relationship between the inputs and the output through a set of parameters: the weights and the bias. By choosing the bias and weights appropriately, the model can represent elementary relationships between data.

2. Multi-Layer Networks: Relations between data are often more complex. The perceptron can be interpreted as one layer of a neural network, sometimes called a single-layer perceptron. An essential concept of deep learning is to build a large (or deep) network of perceptrons arranged in multiple layers, which leads to the multi-layer perceptron model.

Each hidden/output node is now related to all the nodes of the previous layer: we say that the network is fully connected. For each layer, we can write the output/input relation between the layer l + 1 and the layer l:

We have used vector notations in bold to encode the values of x, h, and b at the l-th layer, and denoted W the matrix of weights. Following the figure we have for example h¹ = (h¹₁, h¹₂, h¹₃, h¹₄). The weighted sum from the perceptron is extended to be a matrix-vector product for the multi-layer perceptron, for example

In short, we can think about the deep structure of the multi-layer perceptron as a multiple composition of perceptrons. This gives more unknown weights and bias parameters to be found, but more freedom for the model to adapt to complex data relationships. The action of obtaining the output y from the input x is called a forward pass.

Implementing the forward pass in Python

We can implement the forward pass of a multi-layer perceptron in a few lines. We use random weights and bias for the example:

import numpy as np

sizes = [3, 4, 3, 1]  # number of neurons in each layer
num_layers = len(sizes)  # one input layer, 2 hidden layers and one output layer

# the first layer does not have bias
biases = [np.random.randn(sizes[i], 1) for i in range(1, num_layers)]
# we use the transpose of the weight matrices
weights = [np.random.randn(sizes[i], sizes[i-1]) for i in range(1, num_layers)]

Once we have initialized the network size and its parameters, we are ready to implement the forward pass. We can then feed the network with any arbitrary input.

# Define an activation function
def g(z):
    """Sigmoid function."""
    return 1.0/(1.0 + np.exp(-z)) 

def network_forward(h): # define output/input relation for a MLP
    """Return the output of the network given input h."""
    for b, w in zip(biases, weights):
        z = np.dot(w, h) + b # mat-vect product + bias
        h = g(z) # apply activation function
    return h

# Example input
input = np.array([[1.],[2.],[3.]]) # a column vector (1,2,3)^T
output = network_forward(input)
print(output)

Comparison between a single and multi-layer perceptron

We compare the capabilities of the single and multi-layer perceptron models on a toy classification example. The input dataset consists of points in the two-dimensional plane arranged in a spiral fashion. The output is converted into class probabilities that represent three colors: red, green, and purple. We use two models

• Single-layer Perceptron: 2 inputs (x-y coordinates) and 3 outputs (3 colors)
• Multi-layer Perceptron: 2 inputs (x-y coordinates), one hidden layer with 50 connected neurons, and 3 outputs (3 colors)

The networks are trained to find/learn the “good” values for the weights and biases. Thanks to its hidden layer, the multi-layer perceptron can draw a more complex boundary decision compared to the single-layer perceptron. We gain classification accuracy, but we need to pay the price of training all the parameters introduced by the hidden layer.

We have seen a simple example of a deep neural net architecture and its modeling capability. The same model structure can be used for recognizing handwritten digits. To complete our introduction to deep learning, we need to shed light on how to train the network, that is how to find the best parameters that fit the data. This will be the topic of the second part of this post.

Stay tuned!

In mathematical optimization, solving Mixed-Integer Nonlinear Programming (MINLP) problems pose unique challenges, especially when dealing with complex systems like cooling production. In this analysis, we delve into the performance of three prominent MINLP solvers — Baron, Bonmin, and Couenne — using a concrete example of a cooling production system. Through the use of the optimization framework Pyomo, I propose in this post to explore their speed, reliability, and suitability for convex and non-convex problem sets.

What is mixed-integer nonlinear programming?

Before we dive into the heart of the topic, let’s define briefly what we mean by MINLP. MINLP seeks to minimize a function that depends on both continuous and discrete variables, making it a versatile tool for optimization. Convex MINLP problems enjoy certain advantages, while non-convex ones introduce complexities such as local optima, challenging the solver’s efficiency and accuracy. Here is how it looks in mathematical terms:

Strategy: Algorithms Employed by MINLP Solvers

MINLP solvers leverage a blend of Linear Programming, Integer Programming, and Nonlinear Programming algorithms. Two classical methods, Branch-and-Bound (BB) and Outer Approximation (OA), are commonly used to address MINLPs.

Branch-and-Bound (BB): A foundational technique that iteratively partitions the problem space into smaller sub-problems, converging towards an optimal solution.

In the illustration, the partition of the problem into sub-problems is represented with a tree structure where each of these new constraints represents a branching decision. The tree is terminated until the solution of the continuous relaxation of the sub-problem is either (i) infeasible, (ii) integer feasible or (iii) lower bound value is not anymore smaller than the value of the incumbent solution.

Outer Approximation (OA): Utilizes decomposition and relaxation to solve a sequence of nonlinear programming subproblems and mixed-integer linear master programs alternately. This method adds cuts at solutions of NLP subproblems and uses the solutions of NLP subproblems to generate cuts which are added to the new relaxations. In the end, there should be enough cuts so that the solution of the MIP relaxation will be also valid to the nonlinear constraints.

The flow chart below shows us the approach of the outer approximation in solving the MINLP problem. The idea is to iteratively build a sequence y0, y1 … yk from an initial solution guess y0, and update the solution if the MIP is still feasible.

Convex vs. Non-Convex MINLP

Convex MINLP problems benefit from well-defined properties, while non-convex ones present greater computational challenges due to the presence of local optima. Techniques like Convexification, Convex Envelopes, and Spatial Branch-and-Bound are employed to tackle these complexities effectively.

MINLP Solvers: A Closer Look

Major (Local) MINLP Solvers

Bonmin: This open-source solver offers a range of algorithms catering to both convex and non-convex problems. Notably, its Branch-and-Bound algorithm serves as a heuristic for non-convex scenarios.

Global Solvers

Couenne: Developed on top of Bonmin, Couenne specializes in finding global optima for non-convex MINLPs. Based on a spatial branch-and-bound framework, this solver obtains a lower bound through an LP relaxation using the reformulation techniques. The algorithm is enhanced by several bound-tightening procedures, a separator of disjunctive cuts, as well as a recently introduced feasibility pump heuristic, and symmetry handling.

Baron: A commercial solver renowned for its efficiency in solving non-convex optimization problems to global optimality. Its spatial branch-and-bound technique, coupled with domain reduction strategies, makes it a top choice for complex scenarios.

The table below presents the main features of each solver:

Computational Study:

Test Problem Formulation

An example of a cooling production system with 3 chillers A, B and C.

The key decision for this problem is to find the optimal combination of chiller status (on/off) together with the operational parameters for the chillers: water flows and return water temperature in the primary loop.

Cost function:

The cost function to be optimized represents the total thermal energy consumption, here the sum of the 3 thermal chillers

The cost function takes the form :

Sets:

Variables:

Constraints:

Results

In our computational study of a cooling production system, we evaluate the performance of Baron, Bonmin, and Couenne across varying problem sizes, focusing on resolution time and objective values. While all three solvers employ a branch-and-bound approach, their effectiveness varies across convex and non-convex problem sets.

Observations:

For Convex Problems: Bonmin or Baron are preferable options, with Bonmin being open source and Baron offering commercial-grade accuracy.

For Non-Convex Problems: While Bonmin operates as a heuristic, both Couenne and Baron are exact solvers. However, Baron’s performance on resolution time may exhibit instability, especially in specific cases like the chiller optimization use case.

Resolution Time: Despite being global solvers, the stability of resolution time for Baron may raise concerns, indicating potential variability in performance.

In conclusion, the choice of MINLP solver hinges on the problem’s convexity and the balance between accuracy and computational efficiency. While global solvers like Baron and Couenne offer comprehensive solutions, local solvers like Bonmin provide viable alternatives for convex optimization. By understanding the nuances of each solver, practitioners can make informed decisions to optimize their problem-solving workflows effectively.

References

[1] Amini, Alexander, et al. “Spatial uncertainty sampling for end-to-end control.” arXiv preprint arXiv:1805.04829 (2018).

[2] Pietro Belotti. Couenne: a user’s manual. Technical report, Technical report, Lehigh University, 2009.

[3] Pierre Bonami and Jon Lee. Bonmin user’s manual. Numer Math, 4:1–32, 2007.

[4] Nikolaos V Sahinidis. Baron — branch and reduce optimization navigator. User’s Manual Version, 4, 2000.

Objective: The task was clear: for each month, calculate the average number of rainy days based on the CHIRPS dataset. To accomplish this, I needed to pre-process the data using Google Earth Engine. For me, I decided to use Google Earth Engine’s Python API as I’m more familiar with Python, but there are also JavaScript options available for those who prefer that language.

Environment Setup: Before diving into the implementation, let’s set up our environment. We’ll use Google Colab for this project. First, we need to install the required packages and authenticate with Google Earth Engine.

# Environment settings
!pip install earthengine-api --upgrade
!pip install pycrs
!pip install rasterio
 
# Authentication and initialization
import ee
ee.Authenticate()
ee.Initialize(project='ee-linhha53')

Using Google Earth Engine: Now that our environment is set up, let’s proceed with utilizing Google Earth Engine to pre-process the CHIRPS data.

• Objective: The average number of rainy days per month for the area of interest
• Dataset: CHIRPS DAILY The daily rain data
• Tool: Using Google Earth Engine through Python API

STEP 1: Create an ImageCollection with CHIRPS daily precipitation data.

# Define start and end years, dates (1981 - 2023)
start_year = 1981
end_year = 2023
start_date = ee.Date.fromYMD(start_year, 1, 1)
end_date = ee.Date.fromYMD(end_year, 12, 31)
 
# Initialize the library.
chirps = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY')
           .filterDate(start_date, end_date)

STEP 2: Create a binary mask for rainy days.

Rainy day: precipitation > 0 mm/day. threshold=0

def create_rainy_mask(image, threshold=0):
    """Function to create a binary rainy day mask. """
    rainy_mask = image.select(['precipitation']).gt(threshold)
    return image.addBands(rainy_mask.rename('rainy_mask'))
 
# Map the function over the ImageCollection
rainy_mask = chirps.map(create_rainy_mask).select('rainy_mask')

STEP 3: Calculate the sum of rainy days for a given month in a given year

Create years and months sequences:

# Create an ee.List object with a sequence of integer numbers between a start year and an end year
years = ee.List.sequence(start_year, end_year)
# Create an ee.List object with a sequence of integer numbers between 1 and 12 (months)
months = ee.List.sequence(1, 12)
 
def yearly_monthly_sum(year):
    """Function to map the monthly_sum() to every month in the months sequences."""
    def monthly_sum(month):
        """Function to calculate the sum of rainy days for a given month in a given year."""
        w = rainy_mask.filter(ee.Filter.calendarRange(year, year, 'year')) \
                      .filter(ee.Filter.calendarRange(month, month, 'month')) \
                      .sum()
        return w.set('year', year) \
                .set('month', month) \
                .set('system:time_start', ee.Date.fromYMD(year, month, 1))
    return months.map(monthly_sum)
 
# Map the yearly_monthly_sum() function to every year in the years sequence to get an ImageCollection with total rainy days in each month in each year
monthly_rainy_days = ee.ImageCollection.fromImages(
    years.map(yearly_monthly_sum).flatten()
)

STEP 4: Calculate the average rainy days for a given month across all years

def monthly_mean(month):
    """Function to calculate the average rainy days for a given month across all years."""
    # Round to the nearest integer.
    w = monthly_rainy_days.filter(ee.Filter.eq('month', month)).mean().round()
    return w.set('month', month) \
            .set('system:time_start', ee.Date.fromYMD(1, month, 1))
 
# Map the monthly_mean() function to every month in the months sequence and make an ImageCollection from the obtained images
mean_monthly_rainy_days = ee.ImageCollection.fromImages(
    months.map(monthly_mean).flatten()
)

STEP 5: Clip to the area of interest

# Add the area of interest
aoi = geemap.shp_to_ee("/content/drive/MyDrive/data/aoi.shp")
 
def clip_image_collection(image):
    """Function to clip each image to the area of interest shapefile."""
    return image.clip(aoi)
 
# Map the clip_image_collection() function to the mean_monthly_rainy_days collection 
mean_monthly_rainy_days_clipped = mean_monthly_rainy_days.map(clip_image_collection)

Displaying the Results: Now that we have processed the data, let’s visualize the results.

1. Interactive map:

The geemap library can be used to display ee.Image objects on an interactive ipyleaflet map.

2. Visualizing results:

• Export the raster image to Drive in format Tif File
• Visualizing Tif File Using Matplotlib and GDAL

Conclusion: In conclusion, this project was a rewarding experience that showcased the power and versatility of Google Earth Engine in analyzing geospatial data. By leveraging this platform, I was able to efficiently process the CHIRPS dataset and derive meaningful insights into rainfall patterns for the specified area of interest. Moving forward, I’m excited to explore more projects that harness the potential of geospatial analysis tools for addressing real-world challenges.

Additional Resources:

• Github Codes
• Google Earth Engine
• CHIRPS Dataset

Feel free to reach out if you have any questions or insights to share about this project!