Model with 100+ parameters, is there any use of SBI?

[mateuszlickindorf]

Hello!

I am an undergraduate student and a true beginner in the world of SBI. I have found the SBI library interesting and am currently trying to optimise a model of Auditory Cortex using it. There have been attempts at doing it using an evolutionary algorithm [1], but i wanted to try something different.

The parameters I am trying to optimise are weight matrices describing the connections between certain regions of the auditory cortex. There are two matrices both of size 14x14 but reducing the connections to anatomically reasonable ones i am left with 109 parameters total. The 109 (92 from first and 17 from the second matrix) parameters after running through my (external, MATLAB based) simulator result in a size 251 vector output. I have access to some computational resources, so performing 1M+ of simulations won't be a problem. EDIT: For now i have trained my density estimator at 10K observations which is as far as i understand a very low number.

After doing some reading, I am pretty sure using SNPE for models with 30+ parameters is computationally too heavy to be feasible. How about the other methods? What would be an optimal approach to a multidimensional problem like that? The head of my lab has no experience using SBI or other inference-based approaches, so he can't really help me here. Are there any ways of approaching the multidimensionality problem, maybe by using the other SBI methods that allow sampling the posterior in a different way (more below)?

I am also surprised to find out, that training the neural density estimator went relatively quick, but the sampling from posterior is taking ages. Is it backed by theory, or is it possibly caused by a mistake in code ie. wrong definition of sample dimensions? I would assume sampling from a built posterior although the size of 109 dimensions isn't that computationally heavy, opposing to the training process of a neural network?

The authors of the original auditory cortex model have suggested weight matrix values that when input to the simulator result in a mediocre fit, maybe there is a way to sample the parameter space in the regions relatively close to those suboptimal suggested ones?

Sorry if my questions sound very basic, I did reading of both the documentation as well as additional resources related to this library and am still puzzled by the questions stated above.

Kind regards
Mateusz

Ewelina Tomana, Nina Härtwich, Adam Rozmarynowski, Reinhard König, Patrick J.C. May, Cezary Sielużycki, Optimising a computational model of human auditory cortex with an evolutionary algorithm, Hearing Research, Volume 439, 2023

[mateuszlickindorf]

here is some code i wrote, it might be a bit chaotic due to me copypasting it from Jupyter, but i believe it might be helpful:
Mind you this is the code for a 1000 simulations i have done before just to test the data processing pipeline.

file_path = 'simulated_results__allSOIs_kFF-2_tau2p006_t-50to200ms1000.csv

# lists to store extracted data
theta_list = []
x_list = []

with open(file_path, 'r') as file:
    lines = file.readlines()

# Process in blocks of 7 lines, thus range is number_of_simulation_data_needed * 7
for i in range(0, 7000, 7):
    try:
        simulation_number = lines[i].strip()  #simulation index
        theta_data = lines[i + 1].strip().split(',')  # 109 matrix weights
        soi_data = lines[i + 5].strip().split(',')  # 251 SOI data, SOI i + n where n = SOI number + 1, here SOI 4

        theta_list.append(list(map(float, theta_data)))
        x_list.append(list(map(float, soi_data)))
        
    except IndexError:
        print("The file may not have complete simulation data for every 7-line block.")
        break

# Conversion to Tensors
theta = torch.tensor(theta_list, dtype=torch.float32)
x = torch.tensor(x_list, dtype=torch.float32)

print("Theta tensor shape:", theta.shape) # #output: Theta tensor shape: torch.Size([1000, 109])
print("X tensor shape:", x.shape) # #output: X tensor shape: torch.Size([1000, 251])

prior = torch.distributions.Uniform(low=torch.zeros(109), high=10 * torch.ones(109)) #prior, the matrix weights are in the [0,10] range

inference = NPE(prior=prior)

inference = inference.append_simulations(theta, x)

density_estimator = inference.train()

#output: Neural network successfully converged after 92 epochs.

posterior = inference.build_posterior(density_estimator)

# Extracting the observation from real data

file_name = 'AidaProcessed.csv'
data = pd.read_csv(file_name, header=None)

soi4 = data.iloc[3].values # I need a specific signal among 5 available, SOI4
soi4_tensor = torch.tensor(soi4, dtype=torch.float32) 
print(soi4_tensor.shape) #output: torch.Size([251])

# observed data x_obs is just the soi4_tensor
theta_posterior = posterior.sample((100,), x=soi4_tensor, show_progress_bars=True)

Is my understanding of the training process corect?

[michaeldeistler]

Hey!

1. I don't think that other methods (NRE/NLE) will be more simulation efficient. For 109 parameters you might need sequential methods, see e.g. Truncated proposals for scalable and hassle-free simulation-based inference here.

2. Yes, you will definitely need many simulations. 1M is a good starting point though.

3. In your code, please use

from sbi.utils import BoxUniform

prior = BoxUniform(torch.zeros(109), 10 * torch.ones(109))

instead of your current prior. This might already fix the slow sampling, but if it does not:

4. I would strongly expect that 1000 simulations is just way too few to generate proper results. In a 109 dimensional space, the density estimator will place almost all its mass outside of the prior in at least one of the dimensions. The posterior automatically rejects samples that are outside of the prior, which is probably why you are seeing such slow sampling. To sample from the raw density estimator (i.e. without rejecting samples outside of the prior bounds, just sample from the density_estimator:

density_estimator = inference.train()
samples = density_estimator.sample((1000,), condition=soi4_tensor)

Hope this helps!
Michael