I’ve had success with a similar problem (Nonlinear ODE which is a function of order 10 parameters) by using surrogate modeling techniques (Multidimensional gaussian process - #52 by NateAM). I build a response surface capturing some output as a function of my incoming parameters. You could do dimensionality reduction on your outputs using PCA and on your inputs using CCA to make the surrogate model more friendly. I use the surrogate modelling toolbox (SMT: Surrogate Modeling Toolbox — SMT 0.8.0 documentation) as a blackbox likelihood function. I like to model the log-likelihood directly (multivariate input, scalar output) but you may have different needs.