mc points normalization and qNIPV acquisition function values

[ramseyissa]

[EDITED - question 3]
Hello all,

I just had a few general questions regarding mc points handling and the qNIPV acquisition function.

Background
I am working on a problem that has 8 features and 6 outputs regrading a material science problem.

This is my code for model initialization:

def initialize_model(train_x, train_obj):
    
    print("Initializing model")
    
    input_transform=Normalize(d=train_x.shape[-1])
    models = []
    for i in (pbar0 := trange(train_obj.shape[-1], leave=True,colour="red",position=0)):
        train_y = train_obj[..., i : i + 1]
        models.append(
            SingleTaskGP(
                train_x, train_y, input_transform=input_transform,
                outcome_transform=Standardize(m=1),
            )
        )
        pbar0.set_description(f"Initializing Model {i+1}/{train_obj.shape[-1]}")
        
    model = ModelListGP(*models)
    mll = SumMarginalLogLikelihood(model.likelihood, model)
    # fit_gpytorch_mll(mll)
    print("Fitting model with multiple attempts...")
    fit_gpytorch_mll(
        mll=mll,
        max_attempts=10,
        pick_best_of_all_attempts=True,
    )
    print("Model Fitted!")
    return mll, model

After I return a fitted model, it is passed to my optimizer.

BATCH_SIZE = 2
NUM_RESTARTS = 10
RAW_SAMPLES = 25

def optimize_qnipv_and_get_observation(model, train_x):
    print("Start Of Optimization Process")
    """Optimizes the qNIPV acquisition function, and returns a new candidate and
    observation."""
    
    seeds = [13]
    experiment_results = {}
    
    mcp_lower_bounds = torch.zeros_like(bounds[0])
    mcp_upper_bounds = torch.ones_like(bounds[1])
    # mcp_lower_bounds = original_bounds[0]  # Lower bounds of actual space -unnormalized
    # mcp_upper_bounds = original_bounds[1]  # Upper bounds of actual space -unnormalized
    mcp_bounds = torch.stack([mcp_lower_bounds, mcp_upper_bounds])
    try:
        mcp = draw_sobol_samples(bounds=mcp_bounds, n=1024, q=1, seed=42).squeeze(-2)
        
        mcp = mcp.to(device=device, dtype=torch.double)
        print(f"MC points shape: {mcp.shape}")
        
        weights = torch.ones(y.shape[-1], device=device, dtype=torch.double)
        posterior_transform = ScalarizedPosteriorTransform(weights=weights)
        
        acq_func  = qNegIntegratedPosteriorVariance(
                model=model,
                mc_points=mcp,
                posterior_transform= posterior_transform,
            )
    
        for i in (pbar1 := trange(len(seeds), leave=True,colour="red")):
            set_seeds(seeds[i])
            
            candidates, acq_value = optimize_acqf(
                acq_function=acq_func,
                bounds=original_bounds,
                q=BATCH_SIZE,
                num_restarts=NUM_RESTARTS,
                raw_samples=RAW_SAMPLES,
                options={"batch_limit": BATCH_SIZE, "maxiter": 200},
                sequential=True,
                inequality_constraints=None,  
                equality_constraints=equality_constraints,  
            )
            
            predictions = model.posterior(candidates).mean
            
            experiment_results[seeds[i]] = {
                'candidates': candidates,
                'acq_value': acq_value,
                'predictions': predictions,
            }
            
            candidate_str = "\n".join(
                    [f"C{j+1}: {c.tolist()}" for j, c in enumerate(candidates)]
                )
            pbar1.set_description(
                    f"Seed {seeds[i]} | "
                    f"Acq Values: {acq_value.tolist()}"
                )
                
            print(f"\nCandidates:\n{candidate_str}")
            print(f"Predictions: {predictions.squeeze().tolist()}")
                
    except Exception as e:
        
        print(f"Error occurred: {str(e)}")
        print(f"MC points shape: {mcp.shape if 'mcp' in locals() else 'Not created'}")
        print(f"Model input dimension: {train_x.shape[-1]}")
        raise
    
    return experiment_results

Here lies some confusion on my end.

1. Are the mc points that are generated suppose to be generated in the unnormalized / original design space, since they are being passed to the constructed fantasy model(see code below)?
2. Another question I had was, can the weights be selected in a more constructive manner on my end (maybe identifying some domain knowledge that can be incorporated into the weights)?
3. One thing I noticed in comparing acquisition values for a batch of candidates that are generated from qLogEHVI and qNIPV is that when say two or three candidates are generated using qLogEHVI the acquisition value for each point is in descending order, meaning the first candidate usually has the highest acquisition value. When I compare that to candidates generated from qNIPV for example, i noticed that it is the opposite, the values tend to be in ascending order. For example:
   'acq_value': tensor([-2991.0874, -2976.0832], dtype=torch.float64)
   This could be a misinterpretation on my end, but I was under the impression that since the values are negated, shouldn't i be selecting the candidate that has the closest value to zero (as it would be the maximum)? So given my acquisition values above, wouldn't the better choice be the second candidate given the above list of acquisition values?
   [EDIT] Thinking a little more about this, i believe this is happening due to sequential being set to True so for qNIPV the variance is actually approaching zero with every newly considered candidate given the fantasy model, where as for qLogEHVI later candidates might only find smaller improvements than initial candidates.
4. Finally, given that my problem is materials science related and that objectives often tend vary together, is this problem more appropriate as a multi-task multi-output problem?

Thank you in Advance, and please let me know if you need me to clarify anything on my end!

source code

def forward(self, X: Tensor) -> Tensor:
        # Construct the fantasy model (we actually do not use the full model,
        # this is just a convenient way of computing fast posterior covariances
        fantasy_model = self.model.fantasize(
            X=X,
            sampler=self.sampler,
        )
...

# evaluate the posterior at the grid points
        with settings.propagate_grads(True):
            posterior = fantasy_model.posterior(
                mc_points, posterior_transform=self.posterior_transform
            )

        neg_variance = posterior.variance.mul(-1.0)