We enumerate over the whole domain to compute the logp of the conditionally dependent variables, not all model variables. If there’s an overlap in the set of conditionally dependent variables of two marginalized variables then you get a quadratic number of evaluations yes.
On the plus side it’s all vectorized when possible and you can throw it at a GPU. So that can help somewhat.
We’re also adding QMC marginalization, with which the user can control the number of eval points although that had marginalization of continuous variables in mind and tends to use an ungodly number of evals as well. Not sure if it can be used for discrete variables and if it allows one to grow more linearly with the number of marginalized variables: Marginalize continuous variable via QMC integration by ricardoV94 · Pull Request #353 · pymc-devs/pymc-experimental · GitHub
And we’re planning to add Laplace integration as well (not sure about the exact term, @theo will know better)
Then we have special cases like the HMM where we dispatch on the well known algorithm.
By the way this started as a special model subclass but we’re refactoring it to be a model transformation you can apply on vanilla models like other transformations: do, observe. So users won’t have to use a different object on their workflow.