# Optimisation in twinLab

In many application areas, often we are not only interested in creating a design or process that just nominally does the job, but is also the best design by some metric. The process of searching for the best design(s) is called *optimisation*. **twinLab** provides the **Recommend Module** on our cloud-based platform to enable users to solve optimisation problems without requiring either deep technical understanding of optimisation methods or local computational resources.

Optimisation problems can vary significantly in complexity, such as minimising the distance required to travel between two locations, maximising the aerodynamic efficiency of a vehicle, or maximising the efficiency of solar panels while minimising costs. As the complexity of the optimisation problem increases, implicit design (by experience, judgment, opinion) becomes more difficult, with increased risk of not obtaining the true optimum. One of the features provided by twinLab is **Bayesian Optimisation**: a powerful algorithmic technique for estimating the global optimum of a function, without making any assumptions about the form of said function. Bayesian Optimisation is a probabilistic method that allows uncertainty quantification, so that the risks associated with a decision can be clearly defined.

The metric with which optimality can be judged is often referred to as the *objective* or *objective function*. Given some continuous objective function, denoted as

$$
\begin{align}
f &: x \rightarrow y, & x \subset \mathbb{R}^n, y \subset \mathbb{R},
\end{align}
$$

The aim of global optimisation is to find either

$$
x^\ast = \arg\max_x f(x)
$$

or

$$
x^\ast = \arg\min_x f(x)
$$

Generally, *global* optimisation is much more difficult than *local* optimisation, because simply determining whether the current best estimate is a local extrema is very challenging.

In the ideal scenario, the objective function has known analytical form, and the optimum can be mathematically determined. However, in many applications in both theoretical and practical domains, the exact analytical form of the objective function is unknown. Here, exploiting what information is available about the unknown objective, such as linearity, convexity, or differentiability, can lead to efficient algorithms. If no information at all is available, then one can try to obtain a very large number of evaluations of the objective over the input space, with the aim of exactly replicating the objective function and then determining the optima.

*Finding a global optimum of a function using sparse, unreliable data*

In contrast to the above scenarios, Bayesian Optimisation is best applied when the structure of the objective function is either completely unknown, unreliable, or impractical to obtain. This may arise due to a number of challenges:

- **High computational cost**: data used to evaluate the objective function are typically obtained from computer simulations, laboratory experiments, or real-world simulations. Here, constraints in time, cost, or materials result in a very small number of possible objective evaluations.
- **Unreliable data**: the evaluations of the objective function may contain irreducible uncertainty stemming from various sources of measurement or instrument errors. Computer simulation results may provide insufficient significant digits, and sensor measurements may be noisy.

Uncertainty arising from noise and lack of data greatly increases the difficulty of solving global optimisation problems.

# Overview

*Bayesian Optimisation loop, consisting of an outer loop and an inner loop*

Bayesian Optimisation provides an elegant framework to solve global optimisation problems with uncertainty. There are two important components:
- **Surrogate model**: First, we create a surrogate model (commonly referred to as the **emulator**) to estimate the objective function based on the currently available (noisy) data, and capture the uncertainty of this estimate.
- **Acquisition function**: We design an *acquisition function* that intelligently selects where to sample next. This is **in itself a global optimisation problem** (referred to as the inner loop). However, optimising the acquisition function is much easier by using the surrogate model as an inexpensive replacement of the expensive objective function.

Each iteration of the Bayesian Optimisation loop involves evaluating the two components above.

## The Surrogate Model

The surrogate model represents an approximation of the objective function based on the currently available data, and typically it is significantly less expensive to evaluate the surrogate than to evaluate the objective function.

One of the most commonly used class of surrogate models that incorporate uncertainty are **Gaussian Processes** (GPs). They are not the focus of the current article, but interested readers may want to refer to **[INSERT OTHER ARTICLE LINK]** for a technical description. Here, we refer to the surrogate simply as

$$
\tilde{f}(x) = \mathcal{GP}(x) = \mu(x), \sigma(x)
$$

GP surrogates operate under the Gaussian assumption, where $\mu(x)$ denotes the model mean and $\sigma(x)$ denotes the model standard deviation evaluated at an input $x$.

Initially, $\tilde{f}(x)$ may be constructed using no data, as objective evaluations are not yet available before the optimisation process.

*Initial surrogate model with no data*

In the absence of information, any sample location is equally likely, and one may opt to obtain an initial dataset

$$
\mathcal{D}_0 = \{x_i, y_i\}_{i=0}^{N_0}
$$

where $N$ denotes some initial sampling budget. This may be obtained using random sampling, or existing data may be available from previous experiments, or a number of other strategies.

*Surrogate model with an initial dataset $\mathcal{D}_0$*

Here, we observe that determining a maximum is non-trivial: the mean estimate of $\tilde{f}(x)$ is already reasonably close to the true maximum, however there is very large uncertainty around $x = 0.4$. Without our knowledge of the (unknown) objective function, it is difficult to decide where to sample next. Here, there are two distinct approaches:

- **Exploitation**: continue to sample around the current best estimate of $x=0.75$, to try to improve further.
- **Exploration**: sample around $x = 0.4$, where uncertainty is currently highest, in order to obtain more information about the objective.

In the Bayesian Optimisation setting, we now select and optimise an *acquisition function*. Many known acquisition functions aim to achieve a balance between exploitation and exploration.

## The Acquisition Function

An acquisition function aims to quantify the expected utility of evaluating the objective function on any particular input $x$ in the search space. In the following, we examine the behaviour of a few of the most common acquisition functions. For simplicity, the notation here assumes the objective is maximising $f(x)$ with Gaussian likelihoods.

### Upper Confidence Bound

$$\text{UCB}(x\,; \beta) = \mu(x) + \beta\sigma(x)$$

The Upper Confidence Bound (UCB) function involves a linear sum of the surrogate model mean and uncertainty at a given input $x$. The hyperparameter $\beta$ the user to decide the balance between exploitation and exploration. Increasing $\beta$ incentivises exploration, whereas decreasing $\beta$ towards 0 incentivises exploitation.

### Probability of Improvement

$$
PI(x;\xi) = \Phi\left(\frac{\mu(x) - y_\max - \xi}{\sigma(x)}\right)
$$

The Probability of Improvement (PI) function yields the probability that evaluating the objective surrogate at a given $x$ will result in an improvement over the current optimum $y_\max = \max_y y_i \in \mathcal{D}$. Here, $\Phi$ denotes the Cumulative Distribution Function of the standard normal distribution. PI is aimed primarily at exploitation, however the hyperparameter $\xi$ allows increased exploration.

### Expected Improvement

$$
EI(x;\xi)=\mathbb{E}[\max\{y_{max}-\tilde{f}(x)-\xi,\, 0\}]
$$

The Expected Improvement (EI) function yields the expected *magnitude* of improving upon the current best estimate $y_\max$, for a given input $x$. For a single sample candidate, the analytical solution can be written as:

$$
EI(x;\xi)= (\mu(x)-y_\max-\xi)\,\Phi\left(\frac{\mu(x) - y_\max - \xi}{\sigma(x)}\right)+\sigma(x)\Phi\left(\frac{\mu(x)-y_\max-\xi}{\sigma(y_\max)}\right)
$$

Here, $\Phi$ denotes the Cumulative Distribution Function of the standard normal distribution. EI achieves a more balanced approach between exploitation and exploration compared to PI, and this can be tuned further with the hyperparameter $\xi$.

### Visualising the Acquisition Functions


The above figure demonstrates (from top to bottom): the current surrogate model $\tilde{f}(x)$, the acquisition functions UCB, PI, and EI evaluated over the (one-dimensional) search space. We observe that the each function assigns different utility values for any particular $x$. Using the exploitation-focused PI function may result in reaching a locally optimal solution faster, whereas the more balanced EI function may result in slower convergence but with less risk of being trapped in a local optimum. For UCB, two evaluations of the acquisition function with different values of the hyperparameter $\beta$ are shown, to illustrate its effect.

The star markers represent the maxima in each acquisition function, which denotes the point in search space providing the maximum utility. Obtaining this location is in itself a global optimisation problem, which is straightforward to perform since the surrogate model should be fast to evaluate. This is then used as the next location to sample the objective function, either via computer simulation or physical measurement. In this manner, the Bayesian Optimisation loop intelligently decides where to sample in each iteration, resulting in a sample-efficient optimisation procedure.

## Putting it all together

Given an unknown objective function

$$
\begin{align}
f &: x \rightarrow y, & x \subset \mathbb{R}^n, y \subset \mathbb{R},
\end{align}
$$

which is expensive or difficult to evaluate, and an initial dataset

$$
\mathcal{D}_0 = \{x_i, y_i\}_{i=0}^{N_0}
$$

we aim to find an $x^\ast$ such that

$$
x^\ast = \arg\max_x f(x)
$$

This involves the outer Bayesian Optimisation loop, with each iteration involving the following steps:

1. Train surrogate model on current dataset $\mathcal{D}$
2. Recommend new sample location(s) by choosing and optimising an acquisition function $\mathcal{A}$. This process is called the inner optimisation loop.
3. Sample the objective function at the recommended point(s) either by computer simulation or physical measurement, and add the new data to $\mathcal{D}$

Within the inner optimisation loop, the acquisition function is denoted as

$$
\begin{align}
\mathcal{A}( x; \tilde{f}; \theta): \mathbb{R}^n &\rightarrow \mathbb{R}
\end{align}
$$

where $\theta$ denotes the vector of hyperparameters used to tune the behaviour of the acquisition function. The inner loop solves the following optimisation problem:

$$
x^* = \arg\max_{x} \mathcal{A}(x; \tilde{f}; \theta)
$$

This optimisation can be performed using standard numerical techniques.

# Code Example

Our `twinlab` Python package provides a very straightforward syntax aimed at allowing users to write a Bayesian Optimisation loop from start to finish with very minimal required code.

# import twinlab
import twinlab as tl

import numpy as np
import pandas as pd
import matplotlib as mpl

from matplotlib  import pyplot as plt

# seed for reproducibility
np.random.seed(12345)

1. Define an (arbitrary) objective function

# true objective
def f(x):
    return -((8*x - 2)**2 *  np.sin(12*x - 4) - 3 * np.sin(5 * np.pi * x) - 2 * (x - 0.5)**2)
    
# Create 100 lineaarly spaced x values between 0 and 1
x = np.linspace(0, 1, 100)

# Evaluate y across the input space
y = f(x)

2. Evaluate initial samples

n_samples = 7

# Select random initial sample locations
x_data = np.random.uniform(0, 1, size=n_samples)

# Evaluate y at sample locations
y_data = f(x_data)

**Note that for simplicity, we have not added noise to the data samples. In practical applications, data always involve some degree of noise.**

# Plot objective function and samples
plt.plot(x, y_true, label="Objective Function")
plt.scatter(x_samples, y_samples, marker='x', color='k', label = "Samples")
plt.show()

*Initial data samples*

3. Define outer loop

def outer_loop_iteration(df):
    """A function defining a single iteration of our Bayesian Optimisation loop.
    
    Parameters
    ----------
    df : pandas.DataFrame
        The current dataset
        
    Returns
    -------
    pandas.DataFrame
        The next dataset 
    """

    # Create twinlab dataset instance and upload our local data
    dataset = tl.Dataset("bayesian_optimisation")
    dataset.upload(df)
    
    # Train an emulator in the dataset df
    emulator = tl.Emulator("emulator")
    emulator.train(dataset = dataset, inputs = ["x"], outputs = ["y"], test_train_split=len(df))
    
    # Run the inner loop to get 1 candidate point 
    x_new = emulator.recommend(1, 'ExpectedImprovement')['x'].to_numpy()

    # Evaluate objective on recommended point(s)
    df_new = pd.DataFrame({'x': x_new, 'y':f(X)})

    # Append to dataset and return result
    df  = pd.concat([df, df_new], ignore_index=True)

    return df

4. Run an iteration of the loop and run inference

# Construct dataset
df = pd.DataFrame({'x': x_data, 'y': y_data})

# Run iteration
df = outer_loop_iteration(df)

# Get inference results
df_test = pd.DataFrame({'x': x, 'y': y})
df_mean, df_stdev = emulator.predict(df_test)

The inference results are shown below:
- The first figure shows the initial surrogate model and the the data sample recommended by the inner optimisation loop, obtained by calling `twinlab.emulator.recommend()` using the Expected Improvement acquisition function (iteration 1)
- In the second figure, a new surrogate model conditioned on the newly recommended data sample is shown alongside the next data sample (iteration 2)
- The third figure shows the surrogate model conditioned on the data points from iterations 1 and 2.

*Iteration 1: initial surrogate model with resulting recommended sample location*

*Iteration 2*

*Surrogate model trained on iteration 2 dataset*

Here, we observe that the Bayesian Optimisation algorithm has reached the true objective maxima in two iterations of the loop. It is pertinent to now discuss the termination conditions of this process. With our knowledge of the objective function, we can safely stop at this point and use the optimum of the surrogate model mean as the solution. However, in the practical setting one may desire to reduce the model uncertainty much further before stopping the optimisation process.

# Conclusion

Bayesian Optimisation is an elegant process to solve difficult global optimisation problems, which excels in drawing uncertainty-aware conclusions when data is not readily available or noisy. The process presented in this article is straightforward, however in practice there are often other factors to consider in order to achieve the best result, including:
- How to handle multi-objective optimisation problems, with mixed maximisation/minimisation criterions
- How to obtain multiple recommended sample locations at once (useful when parallelisation is readily available in simulation or physical measurements)
- How to ensure constraints imposed upon a solution are adhered to
- How to integrate prior knowledge on the noise level in the objective evaluations
- Which optimisation routine to use for the inner loop

Within **twinLab**, most of these factors are by default abstracted away from the user, so that they do not need to be concerned with all of the data science specifics that exist at a deeper level than is presented in this article. Nonetheless, the ability to tune the behaviour of the algorithm to address all of the above factors are still exposed to the user, with clear documentation and examples to guide their specific use cases.