Adversarial Machine Learning

  • “Can my classifier easily be fooled?”
  • Basis: MNIST handwritten digit database.
  • Given: a (small) trained classifier and one well classified image \(\bar x\).
  • Goal: Construct an image \(x\) close to \(\bar x\) that is classified with a different label.

A Multinomial Logistic Regression Classifier

In this example we use multinomial logistic regression to classify digits

clf = LogisticRegression(C=50.0 / 500, penalty="l1", solver="saga", tol=0.1, random_state=4)
pipeline = make_pipeline(clf)
pipeline.fit(X, y)
Pipeline(steps=[('logisticregression',
                 LogisticRegression(C=0.1, penalty='l1', random_state=4,
                                    solver='saga', tol=0.1))])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

A Training Example

  • \(\bar x\) is a grayscale image of \(28\times 28\) (\(=784\)) pixels

  • The output is a vector \(y\) of length 10 (each entry corresponding to a digit)

  • The image is classified according to the largest entry of \(y\)

  • The entries of the vector are computed using SoftMax and are probability estimates


Predicted label is 1

y = [ 0.  98.3  0.1  0.3  0.   0.4  0.1  0.   0.9  0. ]

Multinomial Logistic Regression and SoftMax Function

The multinomial regression computes one regression vector for each label:

\[ r_l(x) = a_l^T x + b_l, \, l=0,\ldots,10 \]

Then the SoftMax function is applied to obtain the vector \(y\):

\[ y_l = g_l(x) = \frac{e^{r_l(x)}}{\sum_{k=0}^{10} e^{r_k(x)}}, \, l=0,\ldots, 10 \]

Tip

\(g\) is typically a function that is now easy to write in Gurobi 12!

A Variant of the Classical Adversarial Problem

  • Index \(l= 1\) of \(y\) is correct. Can we trick the classifier to classify a nearby image with index \(w= 5\)?
  • An image \(x\) with \(y_w > y_l\) is a counter example
  • Here we additionally require the probability of \(l\) to be smaller than \(0.25\)

\[\begin{alignat}{3} \max_{0 \le x \le 1, y} \quad & y_w\\ \mbox{s.t.} \quad & \|x - \bar{x}\|_1 \le \delta\\ & g(x) = y \\ & y_l \le 0.25 \end{alignat}\]

Notations

\(g\) the prediction function of the neural network

\(\delta\) the width of the neighbourhood

The gurobipy Model

m = gp.Model()
delta = 6
right_label = 1
wrong_label = 5

# Create input/output variables
x = m.addMVar(xbar.shape, lb=0.0, ub=1.0, name="x")
y = m.addMVar(10, lb=-np.inf, name="y")

# Constrain x to 1-norm neighbourhood around \bar x
eta = m.addMVar(xbar.shape, lb=0, ub=1, name="abs_diff")
m.addConstr(eta >= x - xbar)
m.addConstr(eta >= -x + xbar)
m.addConstr(eta.sum() <= delta);

# Upper bound on probability of right label
y[right_label].UB = 0.25

# Objective: Maximize the probablilty of wrong label
m.setObjective(y[wrong_label], gp.GRB.MAXIMIZE)

# Insert the neural network in the gurobipy model
pred_constr = gml.add_predictor_constr(m, pipeline, x, y, predict_function="predict_proba")

Solving the Model

Gurobi Optimizer version 12.0.0 build v12.0.0rc0 (mac64[arm] - Darwin 23.6.0 23H222)

CPU model: Apple M3 Pro
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 1589 rows, 1588 columns and 3940 nonzeros
Model fingerprint: 0xa32573da
Model has 20 general nonlinear constraints (20 nonlinear terms)
Variable types: 1588 continuous, 0 integer (0 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [2e-01, 1e+00]
  RHS range        [4e-03, 6e+00]
Presolve model has 20 nlconstr
Added 11 variables to disaggregate expressions.
Presolve removed 1442 rows and 966 columns
Presolve time: 0.01s
Presolved: 248 rows, 634 columns, 4960 nonzeros
Presolved model has 10 bilinear constraint(s)
Presolved model has 10 nonlinear constraint(s)
Warning: Model contains variables with very large bounds participating
         in nonlinear terms.
         Presolve was not able to compute smaller bounds for these variables.
         Consider bounding these variables or reformulating the model.


Solving non-convex MINLP

Variable types: 634 continuous, 0 integer (0 binary)

Root relaxation: objective 1.000000e+00, 282 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    1.00000    0   20          -    1.00000      -     -    0s
     0     0    1.00000    0   16          -    1.00000      -     -    0s
     0     0    1.00000    0   14          -    1.00000      -     -    0s
     0     0    1.00000    0   14          -    1.00000      -     -    0s
     0     0    1.00000    0   14          -    1.00000      -     -    0s
H    0     0                       0.7396035    1.00000  35.2%     -    0s
     0     2    1.00000    0   14    0.73960    1.00000  35.2%     -    0s

Explored 360 nodes (3690 simplex iterations) in 0.20 seconds (0.37 work units)
Thread count was 8 (of 8 available processors)

Solution count 1: 0.739604 

Optimal solution found (tolerance 1.00e-04)
Best objective 7.396035216634e-01, best bound 7.396198121843e-01, gap 0.0022%

Solution

Label probabilites for perturbed image:
 [ 0.  24.4  0.1  0.1  0.  74.   0.   0.   1.4  0. ]

Softmax Model in Gurobi ML

# We want to write y_j = e^z_j / sum_j=1^k e^z_j
# y_j are the output variable z_j = input @ coefs + intercepts

# Store the e^z_j in a nonlinear expression
exponentials = nlfunc.exp(linear_predictor)

# The denominator is the sum over the first axis
denominator = exponentials.sum(axis=1)

# Voila!
gp_model.addConstr(
    output == exponentials / denominator[:, np.newaxis]
)

Gurobi 11 Style Model

linear_predictor_vars = gp_model.addMVar(
    output.shape, lb=-gp.GRB.INFINITY, name="linear_predictor"
)
gp_model.addConstr(linear_predictor_vars == linear_predictor, name="linreg")
predictor_model.linear_predictor = linear_predictor_vars

exponentials = gp_model.addMVar(output.shape)
exponentials = exponentials
denominator = gp_model.addMVar((output.shape[0]), lb=epsilon)

num_gc = gp_model.NumGenConstrs

for index in np.ndindex(output.shape):
    gp_model.addGenConstrExp(
        linear_predictor_vars[index],
        exponentials[index],
        name=predictor_model._indexed_name(index, "exponential"),
    )
gp_model.update()
try:
    attributes = predictor_model.attributes
except AttributeError:
    attributes = {}
for gen_constr in gp_model.getGenConstrs()[num_gc:]:
    for attr, val in attributes.items():
        gen_constr.setAttr(attr, val)
gp_model.addConstr(denominator == exponentials.sum(axis=1))

gp_model.addConstrs(
    output[i, :] * denominator[i] == exponentials[i, :]
    for i in range(output.shape[0])
)

Comparisons of Models Performance

Using Gurobi 12 beta:

  • PWL (aka Gurobi 10 way)
  • Univariate (aka Gurobi 11)
  • Expressions (aka Gurobi 12)
v10 v11 v12
#models # viol. Speedup # viol. Speedup # viol.
50 15 19.5 0 35.2 0

With better perormance try adding one small dense NN layer:

v11 v12
#models # timeout # viol. # timeout # viol. Speedup
29 16 0 1 2 102