> For the complete documentation index, see [llms.txt](https://zavalab.gitbook.io/statschemebook/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://zavalab.gitbook.io/statschemebook/chapter-7-decision-making-under-uncertainty/11-bayesian-optimization.md).

# Bayesian Optimization

* Bayesian Optimization (BO) is a special type of stochastic optimization algorithm.
* BO provides “closed-loop” paradigm that integrates experimental design/data collection, modeling/learning, uncertainty quantification, and decision-making.
* Consider system:

  $$
  Y(u)=\varphi(u)+\epsilon
  $$
* $$\varphi(u)$$ is performance, $$\epsilon$$ is random noise, $$u$$ is decision, and $$Y$$ is observed output.
* Decision $$u$$ represents conditions under which we can run the system (e.g., experimental design).
* Goal of BO is to find conditions $$u$$ that maximize system performance:

  $$
  \max\_{u\in \mathcal{U}};; \varphi(u)
  $$
* Challenge is that we *do not know* system function $$\varphi(u)$$ (system is *black box*).
* We need to build a *data-driven* model using data $$(u\_\omega,y\_\omega),; \omega\in \mathcal{S}$$.
* This model is known as a *surrogate* model and is used to predict performance.
* Can use diverse models (e.g., linear regression, physics, neural nets, basis).
* Popular modeling approach used in BO is *kriging*, as it captures nonlinear behavior and quantifies uncertainty.
* Imagine we have a kriging model built using data $$(u\_\omega,y\_\omega),; \omega\in \mathcal{S}$$.
* Kriging provides mean prediction model $$\mu\_\ell(u)$$ and uncertainty model $$\sigma\_\ell(u)$$.
* Define initial dataset as $$\mathcal{D}*\ell$$ and initial model as $$(\mu*\ell(u),\sigma\_\ell(u))$$.
* Here, $$\ell$$ is an iteration counter.
* We can use prediction/uncertainty model for a couple of tasks:
  * Compute decision $$u\_{\ell+1}$$ that maximizes the expected performance:

    $$
    \mathop{\max}*{u\in \mathcal{U}};\mu*\ell(u).
    $$
  * Compute decision $$u\_{\ell+1}$$ that maximizes information:

    $$
    \mathop{\max}*{u\in \mathcal{U}};\sigma*\ell(u).
    $$
* 1st task is *exploitation*, as it exploits knowledge ($$\mu\_\ell(u)$$) to max performance.
* 2nd task is *exploration*, as it explores input space to max info ($$\sigma\_\ell(u)$$).
* Select new decision $$u'$$ that balances exploration/exploitation by solving:

  $$
  \max\_{u\in \mathcal{U}}; \mu\_\ell(u) + \kappa \sigma\_\ell(u).
  $$
* Objective is known as the *acquisition function* $$AF\_{\ell}(u)=\mu\_\ell(u) + \kappa\cdot \sigma\_\ell(u)$$.
* $$\kappa\in \mathbb{R}\_+$$ is parameter that balances exploration/exploitation.
* Run experiment $$u'$$ in the system to obtain new data point $$(u',y')$$ and augment dataset $$\mathcal{D}*{\ell+1}\leftarrow \mathcal{D}*\ell\cup (u',y')$$.
* Use new dataset to re-train kriging model $$\mu\_{\ell+1}(u),\sigma\_{\ell+1}(u)$$ and use this to obtain new input that maximizes $$AF\_{\ell+1}(u)$$.
* Loop is repeated over iterations $$\ell,\ell+1,\dots$$.
* BO selects experiments to progressively improve model and maximize performance.
* BO captures how humans navigate world by collecting data, extracting knowledge from data (learning models), and leveraging knowledge to make decisions.

## Example: BO for Experimental Design

{% hint style="danger" icon="flask" %}

* A team of scientists want to determine optimal conditions to extract pectin from sunflower heads using oxalic acid as a solvent.
* Need to identify how pectin yield ($$Y$$) is influenced by $$u\_1$$ (extraction temperature), $$u\_2$$ (extraction time), and $$u\_3$$ (oxalic acid concentration).
* Unfortunately, there is *no physical understanding* of system and yield thus needs to be optimized directly using data.
* Apply BO to automatically identify experimental conditions that maximize yield.
  {% endhint %}

{% hint style="danger" icon="flask" %}

* Construct an initial kriging model with 3 experiments.
* Subset of experiments proposed by BO and evolution of yield are shown below.
* BO first runs the 3 inputs at their highest values to reduce uncertainty of model.
* BO then progressively lowers the values of the inputs but eventually realizes that lowering temperature ($$u\_1$$) is not a good strategy.
* As such, it decides to increase temperature while continuing to lower time ($$u\_2$$) and oxalic acid concentration ($$u\_3$$).
* Yield is improved quickly and it eventually saturates, converging to the conditions $$u\_1=1$$, $$u\_2=0.426$$, and $$u\_3=0.24$$.
* BO identifies optimal conditions in less than 30 experiments.
  {% endhint %}

{% hint style="danger" icon="flask" %}
$$
\left\[\begin{array}{ccc}
1.0000 & 1.0000 & 1.0000 \\
1.0000 & 1.0000 & 1.0000 \\
0.9946 & 0.9946 & 0.9692 \\
0.8355 & 0.5647 & 0.5133 \\
0.8697 & 0.7806 & 0.3191 \\
1.0000 & 0.6770 & 0.2967 \\
\vdots & \vdots & \vdots \\
1.0000 & 0.4505 & 0.2219 \\
1.0000 & 0.4426 & 0.2283 \\
1.0000 & 0.4356 & 0.2338 \\
1.0000 & 0.4300 & 0.2380 \\
1.0000 & 0.4259 & 0.2411 \\
1.0000 & 0.4259 & 0.2411
\end{array}\right]
$$

*Subset of experiment history suggested by BO.*
{% endhint %}

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-3b1a1f084c89a94db224a27120bc4257f67072c7%2Fch7_bayesopt_extraction.png?alt=media" alt="" width="420"><figcaption><p>Evolution of yield using experiments proposed by BO.</p></figcaption></figure>

#### Python implementation

{% hint style="success" %}
Notebook [`ch7_bayesopt_yield.ipynb`](https://colab.research.google.com/github/vzavala/StatsBook/blob/master/Python/ch7_bayesopt_yield.ipynb) — open it in Colab and run the cells alongside this section.

MATLAB live script: [`ch7_bayesopt_yield.mlx`](https://github.com/vzavala/StatsBook/blob/master/Matlab-Live/ch7_bayesopt_yield.mlx)
{% endhint %}

Bayesian optimization of pectin extraction (Victor Z and Leo Gonzalez, UW-Madison, 2024).

| Slides                                           | Python                           | Meaning                                         |
| ------------------------------------------------ | -------------------------------- | ----------------------------------------------- |
| $$u=(u\_1,u\_2,u\_3)\in\[-1,1]^3$$               | `u`                              | scaled temperature, time, acid concentration    |
| $$\varphi(u)$$                                   | `phi(u)`                         | yield (a quadratic model plays the "black box") |
| $$\mathcal{D}\_\ell$$                            | `data_u`, `data_y`               | experiments so far                              |
| $$m\_\ell(u)$$, $$\sigma\_\ell(u)$$              | `gp.predict(u, return_std=True)` | kriging mean and SD                             |
| $$AF\_\ell(u)=m\_\ell(u)+\kappa\sigma\_\ell(u)$$ | `AF(u)`, `kappa`                 | acquisition function ($$\kappa=0.2$$)           |

**Note:** This notebook reproduces Example 7.11 (Section 7.7) of *Statistics for Chemical Engineers*: Bayesian optimization (BO) of an extraction system, using a GP model and an acquisition function to choose each new experiment.

```python
import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF
from scipy.optimize import minimize, Bounds

b = np.array([16.5780, 3.8563, 0.9250, 0.4212, -1.2840, -0.8515, -1.0940, -0.1425, 0.1900, -0.1675])

def phi(u):
    u = np.atleast_2d(u)
    u1, u2, u3 = u[:, 0], u[:, 1], u[:, 2]
    return (b[0] + b[1]*u1 + b[2]*u2 + b[3]*u3 + b[4]*u1**2 + b[5]*u2**2 + b[6]*u3**2
            + b[7]*u1*u2 + b[8]*u1*u3 + b[9]*u2*u3)
```

Start from 3 experiments and iterate: fit the kriging model, maximize $$AF\_\ell(u)$$, run the new experiment, and add it to the data set.

```python
kappa, n_iter = 0.2, 20
bounds = Bounds(-np.ones(3), np.ones(3))
data_u = np.array([[0, 0, 0], [1, 1, 1], [-1, -1, -1]], dtype=float)
data_y = phi(data_u)
# alpha adds a small jitter to the kernel matrix: the yields have no noise, and without it the fit is ill-conditioned
gp = GaussianProcessRegressor(kernel=RBF(length_scale=np.ones(3)), normalize_y=True, alpha=1e-4)

for ell in range(n_iter):
    gp.fit(data_u, data_y)
    def neg_AF(u):
        m, s = gp.predict(np.atleast_2d(u), return_std=True)
        return -(m[0] + kappa * s[0])
    u_new = minimize(neg_AF, data_u[np.argmax(data_y)], method='SLSQP', bounds=bounds).x
    data_u = np.vstack([data_u, u_new])
    data_y = np.append(data_y, phi(u_new))

def descale(u, lb=np.array([70, 1, 0.1]), ub=np.array([100, 2, 0.3])):
    """Map [-1, 1] to temperature [C], time [h], concentration [M]."""
    return lb + (u + 1) / 2 * (ub - lb)

best = np.argmax(data_y)
print("best scaled conditions u =", np.round(data_u[best], 3), "  yield =", round(float(data_y[best]), 3))
print("in physical units:", np.round(descale(data_u[best]), 3))

ell = np.arange(len(data_y))
plt.step(ell, np.maximum.accumulate(data_y), 'k-', where='post', label='best so far')
plt.scatter(ell, data_y, color='w', edgecolor='k', label=r'$y_\ell$')
plt.xlabel(r'Experiment $\ell$'); plt.ylabel(r'Observed yield $y_\ell$')
plt.legend(); plt.grid(True)
plt.show()
```

```
best scaled conditions u = [1.    0.463 0.249]   yield = 19.395
in physical units: [100.      1.732   0.225]
```

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-3cd5939bcd306f240cd23c0862c3c92744627d7a%2Fch7_bayesopt_yield-05-1.png?alt=media" alt="" width="520"><figcaption></figcaption></figure>

**Note:** As in the book, BO balances exploration and exploitation, and the experiment history converges to the best extraction conditions in a few experiments.