[Question] What if the return mean & sem of status_quo for each iteration is different?

[FriendLey]

In the online A/B test. We use the BatchTrial mode and set a status_quo(control group). Because of fluctuations in the online metrics, we usually have different mean and sem of status_quo in each iteration.

Suppose we need to optimize a range parameter $a \in [0.0, 1.0]$. The objective metric is video playback time per capita which can only be obtained through online experiments. We started an online experiment with three groups and the control group's parameter is $a = 0.0$. The experiment went through two iterations, one from Monday to Wednesday and the other from Thursday to Saturday. we get the following data frame:

arm_name	trial_index	a	mean(s)	sem	time
0_0	0	0.3	1600	0	Mon. - Wed.
0_1	0	0.5	1700	0	Mon. - Wed.
Control	0	0.0	1650	0	Mon. - Wed.
1_0	1	1.0	2500	0	Thu. - Sat.
1_1	1	0.8	2600	0	Thu. - Sat.
Control	1	0.0	2300	0	Thu. - Sat.

(For convenience, assuming that sem of each group is 0)

if I directly use Models.BOTORCH(experiment=exp, data=exp.eval()) to generate a model, It will prompt Warning: [WARNING 03-16 14:37:02] ModelBridge: Status quo status_quo found in data with multiple features. Use status_quo_features to specify which to use. But It can be seen from the results that the mean between multiple iterations cannot be directly compared. So I can't specify which to use.

My method is to first calculate the relative difference of each group relative to the control group of the corresponding iteration. After processing, I have the following data:

arm_name	trial_index	a	mean(s)	sem	time
0_0	0	0.3	(1600 - 1650) / 1650 = -0.030	0	Mon. - Wed.
0_1	0	0.5	(1700 - 1650) / 1650 = 0.030	0	Mon. - Wed.
1_0	1	1.0	(2500 - 2300) / 2300 = 0.087	0	Thu. - Sat.
1_1	1	0.8	(2600 - 2300) / 2300 = 0.130	0	Thu. - Sat.
Control	1	0.0	0	0	Thu. - Sat.

Then pass this data to Models.BOTORCH. My question are:

1. How does Ax handle this scenario?
2. Is My method correct? If correct, does Ax consider supporting this method internally?

[Balandat]

Your approach is a good one, relativizing the metrics w.r.t. to the status quo is generally the first step for dealing with non-stationary data (assuming that the relative effect remains reasonably stable over time). You'll want to make sure to properly handle the sems, and potentially also do a bias correction, see e.g. section 2.2. in https://arxiv.org/pdf/1803.06336.pdf.

1. Internally, we handles this scenario in much the same way. If there is reason to believe that the relative effect varies over time, we in addition to relativization use a multi-task model that considers the previous evaluations as "similar but not the same" and learns the correlation between these and the "target task" (which is typically the most recent online experiment).

2. Yes, by and large your method is correct, just be careful about properly relativizing the errors, and check for stationarity of the relativized effect. We can definitely support this method - we have a number of internal helpers that we use for handling this, it's probably time to move these into open-source (cc @ldworkin)

[ldworkin]

Yes! Agreed, we should support this, or at least expose the relevant helper functions and make it easier for the user. I'll make a task.

[FriendLey]

Thanks very much for your response. First, I agreed completely using delta method with bias correction to handle the means and sems. I have some new questions about multi-task procedure of Ax.

As @Balandat said:

If there is reason to believe that the relative effect varies over time, we in addition to relativization use a multi-task model that considers the previous evaluations as "similar but not the same" and learns the correlation between these and the "target task" (which is typically the most recent online experiment).

And also As http://www.jmlr.org/papers/volume20/18-225/18-225.pdf modeling different simulator batches with separate tasks.

My understanding is treating different batches as different tasks. I try to simulate this as the following codes(changing from multi_task). I removed the offline task and treat every batch trial of online as a task. I try to simulate but failed. my code is as follows:

code

from copy import deepcopy
import numpy as np
import pandas as pd
from scipy.stats import norm
import sys
import time

from ax.core.data import Data
from ax.core.experiment import Experiment
from ax.core.observation import ObservationFeatures, observations_from_data
from ax.core.optimization_config import OptimizationConfig
from ax.core.search_space import SearchSpace
from ax.core.objective import Objective
from ax.runners.synthetic import SyntheticRunner
from ax.modelbridge.random import RandomModelBridge
from ax.core.outcome_constraint import OutcomeConstraint
from ax.core.types import ComparisonOp
from ax.core.parameter import RangeParameter, ParameterType
from ax.core.multi_type_experiment import MultiTypeExperiment
from ax.metrics.hartmann6 import Hartmann6Metric
from ax.metrics.l2norm import L2NormMetric
from ax.modelbridge.factory import get_sobol, get_GPEI, get_MTGP
from ax.core.generator_run import GeneratorRun
from ax.plot.diagnostic import interact_batch_comparison
from ax.plot.trace import optimization_trace_all_methods
from ax.utils.notebook.plotting import init_notebook_plotting, render

init_notebook_plotting()

def get_experiment(noise_sd=0.1):
    # 1. Create simple search space for [0,1]^d, d=6
    param_names = [f"x{i}" for i in range(6)]
    parameters=[
        RangeParameter(
            name=param_names[i], parameter_type=ParameterType.FLOAT, lower=0.0, upper=1.0
        )
        for i in range(6)
    ]
    search_space=SearchSpace(parameters=parameters)
    # 2. Specify optimization config
    online_objective = Hartmann6Metric("objective_0", param_names=param_names, noise_sd=noise_sd)
    online_constraint = L2NormMetric("constraint_0", param_names=param_names, noise_sd=noise_sd)
    opt_config = OptimizationConfig(
        objective=Objective(online_objective, minimize=True),
        outcome_constraints=[OutcomeConstraint(online_constraint, op=ComparisonOp.LEQ, bound=1.25, relative=False)]
    )
    # 3. Init experiment
    exp = MultiTypeExperiment(
        name="mt_exp",
        search_space=search_space,
        default_trial_type="trial_0",
        default_runner=SyntheticRunner(),
        optimization_config=opt_config,
    )
    return exp

def run_mtgp(noise_sd = 0.1):
    t1 = time.time()
    ### Do BO with online only
    ## Quasi-random initialization
    exp = get_experiment(noise_sd=noise_sd)
    m = get_sobol(exp.search_space, scramble=False)
    gr = m.gen(n=5)
    exp.new_batch_trial(trial_type="trial_0", generator_run=gr).run()
    ## Do BO
    for b in range(3):
        print('bo3 batch', b, time.time() - t1)
        # Fit the GP
        if b > 0:
            m = get_MTGP(
                experiment=exp,
                data=exp.fetch_data(),
                search_space=exp.search_space,
            )
            gr = m.gen(
                n=5,
                optimization_config=exp.optimization_config,
                fixed_features=ObservationFeatures(parameters={'trial_type': 'trial_'+str(b)}),
            )
        else:
            m = get_GPEI(
                experiment=exp,
                data=exp.fetch_data(),
                search_space=exp.search_space,
            )
            gr = m.gen(
                n=5,
                optimization_config=exp.optimization_config,
            )
        # Generate the new batch
        param_names = [f"x{i}" for i in range(6)]
        exp.add_trial_type("trial_"+str(b+1), SyntheticRunner())
        objective = Hartmann6Metric("objective_"+str(b+1), param_names=param_names, noise_sd=noise_sd)
        constraint = L2NormMetric("constraint_"+str(b+1), param_names=param_names, noise_sd=noise_sd)
        exp.add_tracking_metric(metric=objective, trial_type="trial_"+str(b+1), canonical_name="objective")
        exp.add_tracking_metric(metric=constraint, trial_type="trial_"+str(b+1), canonical_name="constraint")
        exp.new_batch_trial(trial_type="trial_"+str(b+1), generator_run=gr).run()
    ## Extract true objective and constraint at each iteration
    return exp

exp = run_mtgp()

Run above code, Ireceive the following bug:

Traceback (most recent call last):
  File "multi_task.py", line 106, in <module>
    exp = run_mtgp()
  File "multi_task.py", line 76, in run_mtgp
    search_space=exp.search_space,
  File "~/anaconda3/envs/r-environment/lib/python3.7/site-packages/ax/modelbridge/factory.py", line 264, in get_MTGP
    status_quo_features=status_quo_features,
  File "~/anaconda3/envs/r-environment/lib/python3.7/site-packages/ax/modelbridge/torch.py", line 75, in __init__
    optimization_config=optimization_config,
  File "~/anaconda3/envs/r-environment/lib/python3.7/site-packages/ax/modelbridge/base.py", line 149, in __init__
    observation_data=obs_data,
  File "~/anaconda3/envs/r-environment/lib/python3.7/site-packages/ax/modelbridge/torch.py", line 89, in _fit
    observation_data=observation_data,
  File "~/anaconda3/envs/r-environment/lib/python3.7/site-packages/ax/modelbridge/array.py", line 83, in _fit
.........
  File "~/anaconda3/envs/r-environment/lib/python3.7/site-packages/gpytorch/lazy/lazy_tensor.py", line 882, in evaluate
    res = self.matmul(eye)
  File "~/anaconda3/envs/r-environment/lib/python3.7/site-packages/gpytorch/lazy/interpolated_lazy_tensor.py", line 402, in matmul
    right_interp_res = left_t_interp(self.right_interp_indices, self.right_interp_values, tensor, base_size)
  File "~/anaconda3/envs/r-environment/lib/python3.7/site-packages/gpytorch/utils/interpolation.py", line 230, in left_t_interp
    summing_matrix = cls(summing_matrix_indices, summing_matrix_values, size)
RuntimeError: size is inconsistent with indices: for dim 1, size is 1 but found index 1

code around the error.

 size = torch.Size((batch_size, output_dim, num_data * num_interp))
type_name = summing_matrix_values.type().split(".")[-1]  # e.g. FloatTensor
if interp_values.is_cuda:
    cls = getattr(torch.cuda.sparse, type_name)
else:
    cls = getattr(torch.sparse, type_name)
summing_matrix = cls(summing_matrix_indices, summing_matrix_values, size)

For the process of modeling each trial as a task, Is This correct? If not, What should I do?

[Balandat]

Hmm I believe this is because the data you get from the experiment does not have any from trials other than trial_0. So the MTGP you're fitting is degenerate in that it's for exactly a single task (task_0), but you're asking it to generate candidates for a new (trial_1) task. You'll need to do sth like exp.new_batch_trial(trial_type="trial_k", generator_run=gr_k).run() for k=1, 2. Note though that in this simple setting the synthetic data will be generated from the same functions, so in this case the best thing would be to just pool the data. If you want to see the MTGP do interesting things, you'll have to manually induce some kind of drift in the metric observations between different k.

[Balandat]

We should probably add some validation to the MTGP model so that we only allow generation for target tasks that are part of the model. cc @qingfeng10

[FriendLey]

Thank you for the explanation. @Balandat

Letting canonical_name of tracking metric equal to name of optimization metric, this error will be fixed.

[bletham]

A general comment:

It's pretty common to see shifts in the metrics like that from one time period to the other, which does introduce a challenge when trying to do multiple iterations of BO using A/B tests.

There are two approaches that we've used to handle this. The first and the easiest is to try relativizing the data for each batch with respect to its control, and see if things are stable across time that way. For example, suppose in my first batch I run (Control, 0_0, 0_1, 0_2, 0_3). Then in the second batch, we will often repeat a few of those arms, and add in a few new ones. So the second batch might be (Control, 0_0, 0_2, 1_1, 1_2) where 1_1 and 1_2 are the new arms, and 0_0 and 0_2 are repeated. (We tend to do large batches, of 20-50 arms, in which case we leave around 5 to be repeated from the initial batch). Then, we compare the "0_0 / control" in the first batch to "0_0 / control" in the second batch. Even though control and 0_0 have both had mean shift, we find that most often the % change from control to 0_0 is relatively stable. If they are stable, then it is safe to compute the % change for all of the arms, and fit a model directly to the % changes across all batches.

There is a diagnostic for checking stability across batches. If you run this command:

from ax.utils.notebook.plotting import render
from ax.plot.diagnostic import interact_batch_comparison

render(interact_batch_comparison(observations=m.get_training_data(), experiment=exp, batch_x=0, batch_y=1))

it will generate a plot that compares the raw values for any repeated arms in batch 0 and batch 1 (as specified by batch_x and batch_y). For instance it might look like this if you have two arms that are present in both batches:

In this case it is good, because we see that the value of the arms were consistent from batch 0 to batch 1 (in this particular experiment this is after relativizing each batch wrt its control). If the points do not line up on the diagonal, then you will likely not get good model fit when combining across batches. It is really important to have a small set (ideally 4 or 5) of arms that are repeated in every batch to be able to verify this stability.

The second approach is to use the multi-task model. With this approach you don't have to relativize the data in each batch since it will effectively learn the adjustment across batches for you. The modeling aspects become a little bit more complicated since, as you discovered, you will now need to specify which batch to generate points for, or plots, which is done with the fixed_features input. Another complication is that if you have many batches, this model will get really slow because it is modeling a separate task for each batch; you'll want to specify a reduced rank. Another thing to note is that when using the multi-task model it is very important (even more important than when relativizing) to have a few arms (e.g. 4 or 5) that are repeated in every batch. This gives the model a fixed reference for learning the adjustment across batches.

Our most typical practice for combining batches in online experiments is the relativization approach, along with repeated arms and verifying with the plot that things are stable after relativization. If things are not stable after relativization, or for other settings like the offline-online, that is when we use the multi-task model.

[FriendLey]

@bletham Thanks for your prompt response.

I have some other questions about Bayesian Optimization in the online A/B test.

If the number of groups is large, then the sample of each group will be less and the variance of metric will be large. In order to reduce the variance of the metric, we could increase the time each batch runs or increase traffic. But time and traffic are the very critical resources of online experiments. My questions are:

1. How to balance the number of groups and traffic?
2. How to decide when to complete the current batch and iterate to the next batch?
3. Are there any other good methods to reduce the variance of a metric？

[bletham]

These are great questions.

1. Basically the question is given a fixed sample budget, do we want more arms with a higher variance per arm, or fewer arms with tighter uncertainty per arm? We have done some simulations on this question which we should probably write up at some point, but I'll give you my opinion based on practical experience. You generally want more arms (with higher variance) than fewer arms. If you have good model fit, then the model will be able to borrow strength across arms, and so the modeled values for those points will have much lower uncertainty than the measured variance. However, it is also important that at least some of the arms that you launch be statistically significantly different. That is, the empirical variance in objective should be low enough that the lowest objective and highest objective arms do not have overlap in their CIs. If the CIs for all of the arms overlap, then "no effect" becomes a valid explanation for the data, and my experience is that the GP will often predict no effect. So, figure out how many samples are required so that the smallest effect and largest effect are significantly different from eachother, and then launch as many arms as possible that all have that many samples.

2. I think this is related to 1. Basically once you have low enough uncertainty to get good model fit, you can optimize and move to the next batch. Another consideration is the stationarity of the outcome and integrating over population shifts / time effects, which is usually our main consideration when deciding when to complete the batch.

3. There are a lot of methods for treatment effect estimation that can reduce variance relative to just taking sample means. If you do some searching for treatment effect estimation or look in a book on experimentation I think you should be able to find lots of information, but the main thing I would point you to that is not terribly difficult to implement is regression adjustment.

[FriendLey]

@ldworkin

Yes! Agreed, we should support this, or at least expose the relevant helper functions and make it easier for the user. I'll make a task.

Any exact timeline on this? We are looking forward to it... Thank you.

[ldworkin]

Hi @FriendLey ! Thanks for the reminder. Sorry for the delay -- we'll bump this up in priority and try to check this code in ASAP.

In the meantime, let me just share a gist with the method I'm planning on adding, so you can let me know if this is what you need: https://gist.github.com/ldworkin/80d80cb8e48001b51b9a61a6ef4810a2

Basically, you should be able to just pass in your dataframe as-is to this, and it should relativize the data for each trial to its respective status quo arm, just like you were doing manually above.

[ldworkin]

This is now in OSS as per above commit, so going to close out. Let us know if you have any other questions!