Global and Preference-based Optimization with Mixed Variables using Piecewise Affine Surrogates (PWAS/PWASp)

Package description
Installation
Basic usage
Contributors
Citing PWAS
License

Package description

We propose a novel surrogate-based global optimization algorithm, called PWAS, based on constructing a piecewise affine surrogate of the objective function over feasible samples. We introduce two types of exploration functions to efficiently search the feasible domain via mixed integer linear programming (MILP) solvers. We also provide a preference-based version of the algorithm, called PWASp, which can be used when only pairwise comparisons between samples can be acquired while the objective function remains unquantified. For more details on the method, please read our paper Global and Preference-based Optimization with Mixed Variables using Piecewise Affine Surrogates.

[1] M. Zhu and A. Bemporad, "Global and preference-based optimization with mixed variables using piecewise aﬃne surrogates," Journal of Optimization Theory and Applications, vol. 204, no. 26, 2025. [bib entry]

Installation

pip install pwasopt

Dependencies:

python >=3.7
numpy >=1.24.3, <2.0.0
scipy >=1.11.1
pulp ==2.8.0
scikit-learn ==1.5.0
threadpoolctl ==3.1.0 (for KMeans from scikit-learn to run properly)
pyparc >=2.0.4
pyDOE ==0.3.8
pycddlib >=2.1.7, <3.0.0
cvxopt ==1.2.7

External dependencies:

MILP solver:

PWAS/PWASp use GUROBI as the default solver to solve the MILP problem of acquisition optimization, which is found to be the most robust during benchmark testing. Alternatively, we also include GLPK, which may introduce errors occasionally depending on the test problem and initial samples. User can also switch to another MILP solver by editing the relevant codes in acquisition.py and sample.py. Check the compatability of the MILP solver with pulp (the LP modeler) at the project webpage.
GUROBI: academic licenses
- configure GUROBI with python
- step-by-step explanation for the installation of Gurobi
GLPK: project webpage
- step-by-step explanation for the installation on Stack Overflow

Basic usage

Examples

Examples of benchmark testing using PWAS/PWASp can be found in the examples folder:

mixed_variable_benchmarks.py: benchmark testing on constrained/unconstrained mixed-variable problems
- Test results are reported in the paper
- Note: to test benchmark NAS-CIFAR10
  - download the data from its GitHub repository
  - indicate the data_path in mixed_variable_benchmarks.py
  - since the dataset is compiled with TensorFlow version 1.x, python version < 3.8 is required (with TensorFlow < 2.x)
other_benchmarks.py: various NLP, MIP, INLP, MIP Benchmarks tested with PWAS/PWASp
- Test results are reported in test_results_on_other_benchmarks.pdf under the examples folder

Case studies

Experimental design with PWAS: ExpDesign

Optimization of reaction conditions for Suzuki–Miyaura cross-coupling (fully categorical)
Optimization of crossed-barrel design to augment mechanical toughness (mixed-integer)
Solvent design for enhanced Menschutkin reaction rate (mixed-integer and categorical with linear constraints)

Illustrative example

Here, we show a detailed example using PWAS/PWASp to optimize the parameters of the xgboost algorithm for MNIST classification task.

Problem discription

Objective: Maximize the classification accuracy on test data. Note PWAS/PWASp optimizes the problem using minimization, and therefore we minimize the negative of classification accuracy.

Optimization variables: $n_c = 4$ (number of continuous variables), $n_{\rm int} = 1$ (number of integer variables, ordinal), and $n_d = 3$ (number of categorical variables, non-ordinal) with $n_{i} = 2$, for $i = 1, 2, 3$. Each categorical variable ($n_{di}$) can be either 0 or 1. The bounds are $\ell_x = [10^{-6}\ 10^{-6}\ 0.001\ 10^{-6}]'$, $u_x = [1\ 10\ 1\ 5]'$; $\ell_y = 1$, $u_y = 10$.

Notes: The 0.7/0.3 stratified train/test split ratio is applied. The xgboost package is used on MNIST classification. The optimization variables in this problem are the parameters of the xgboost algorithm. Specifically, the continuous variables $x_1$, $x_2$, $x_3$, and $x_4$ refer to the following parameters in xgboost, respectively: learning_rate, min_split_loss, subsample , and reg_lambda. The integer variable $y$ stands for the max_depth. As for the categorical variables, $n_{d1}$ indicates the booster type in xgboost where $n_{d1} = {0, 1}$ corresponding to {gbtree, dart}. $n_{d2}$ represents the grow_policy, where $n_{d2} = {0, 1}$ corresponding to {depthwise, lossguide}. $n_{d3}$ refers to the objective, where $n_{d3} = {0, 1}$ corresponding to {multi:softmax, multi:softprob}.

Problem specification in Python

import xgboost
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split
from sklearn import metrics
import numpy as np

# info of optimization variables 
nc = 4  # number of continous variables
nint = 1 # number of integer variables, ordinal
nd = 3  # number of categorical variables, non-ordinal
X_d = [2, 2, 2]  # possible number of classes for each categorical variables

lb_cont_int = np.array([1e-6, 1e-6, 0.001, 1e-6, 1])  # lower bounds for continuous and integer variables
ub_cont_int = np.array([1, 10, 0.99999, 5, 10])  # upper bounds for continuous and integer variables
lb_binary = np.zeros((nd))  # lower bounds for categorical variables, note the dimension is the same as nd, it will be updated within the code
ub_binary = np.array([1, 1, 1]) # upper bounds for categorical variables, note it is (the number of classes-1) (since in the one-hot encoder, the counter started at 0)
lb = np.hstack((lb_cont_int, lb_binary)) # combined lower and upper bounds for the optimization variables
ub = np.hstack((ub_cont_int, ub_binary))

# load dataset
# example code: https://github.com/imrekovacs/XGBoost/blob/master/XGBoost%20MNIST%20digits%20classification.ipynb
mnist = load_digits()  
X, y = mnist.data, mnist.target
X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.7, test_size=0.3, stratify=y,
                                                    random_state=1)  # random_state used for reproducibility
dtrain = xgboost.DMatrix(X_train, label=y_train)
dtest = xgboost.DMatrix(X_test, label=y_test)

# define the objective function, x collects all the optimization variables, ordered as [continuous, integer, categorical]
def fun(x):  
    xc = x[:nc]  # continuous variables
    xint = x[nc:nc + nint]  # integer variables
    xd = x[nc + nint:]  # categorical variables

    if xd[0] == 0:
        mnist_booster = 'gbtree'
    else:
        mnist_booster = 'dart'

    if xd[1] == 0:
        mnist_grow_policy = 'depthwise'
    else:
        mnist_grow_policy = 'lossguide'

    if xd[2] == 0:
        mnist_obj = 'multi:softmax'
    else:
        mnist_obj = 'multi:softprob'
    param = {
        'booster': mnist_booster,
        'grow_policy': mnist_grow_policy,
        'objective': mnist_obj,
        'learning_rate': xc[0],
        'min_split_loss': xc[1],
        'subsample': xc[2],
        'reg_lambda': xc[3],
        'max_depth': int(round(xint[0])),
        'num_class': 10  # the number of classes that exist in this datset
    }

    bstmodel = xgboost.train(param, dtrain)

    y_pred = bstmodel.predict(
        dtest)  # somehow predict gives probability of each class instead of which class it belongs in...

    try:
        acc = metrics.accuracy_score(y_test, y_pred)

    except:
        y_pred = np.argmax(y_pred, axis=1)
        acc = metrics.accuracy_score(y_test, y_pred)

    return -acc  # maximize the accuracy, minimze the -acc

# Specify the number of maximum number of evaluations (including initial sammples) and initial samples
maxevals = 100
n_initil = 20

# default setting for the benchmarks
isLin_eqConstrained = False  # specify whether linear equality constraints are present
isLin_ineqConstrained = False  # specify whether linear inequality constraints are present
Aeq = np.array([])  # linear equality constraints
beq = np.array([])
Aineq = np.array([])  # linear inequality constraints
bineq = np.array([])

Solve use PWAS

One can solve the optimization problem either by explicitly passing the function handle fun to PWAS, or by passing the evaluation of fun step-by step.

Solve by explicitly passing the function handle

from pwasopt.main_pwas import PWAS  

key = 0
np.random.seed(key)  # r

PWAS

Install / Use

README