Good questions even though no single response can address all of your concerns. This will be, at best, a partial answer. It will also be from an applied -- not theoretical level -- it is practical advice. It sounds like you're thinking only of two-way interactions but, obviously, three-way and higher interactions are possible. In fact, one can imagine extending the interactions up to a "fully saturated" limit that would exhaust the degrees of freedom. So, let's limit the discussion to two-way interactions which are easily generalizable to higher ones.
Let's also mention that statistical theory recommends (or prescribes) that the appropriate main effects be retained in any model containing interactions with those effects. This recommendation isn't always followed or observed, but it needs to be noted -- do with it what you will. The more important thing to note is that the interpretation of the model parameters shifts from one of being conditional on the mean values of the features when there are only main effects (i.e., for models that are evaluating the mean as in OLS regression. This would change with differing functional forms such as quantile regression) to an interpretation in the presence of interactions that is conditional on the features being evaluated at zero. In applied terms, this means that, before taking the interaction, the features should be normalized to a mean of zero. Regrettably, none of the points made in this paragraph are without controversy and many statisticians can and will disagree with them. I can provide references supporting these statements if anyone is interested.
That said, the first approach to incorporating interactions can be described as a "global," "brute force" or "shotgun" approach insofar as all possible permutations of (two-way) interactions are evaluated. I say permutations and not combinations since interactions are symmetric -- it doesn't matter if the interaction is F1 with F2 or F2 with F1. Obviously, with more than a few features this approach would quickly result in an unmanageably large number of interactions. Regardless, it is possible and has been done.
A second, somewhat more reasonable approach is to use subject matter expert (SME) knowledge to make judgments as to which interactions are interesting or useful to model. This has the advantage of immediately reducing the number of interactions to a manageable size as well as minimizing time wasted on potentially fruitless exploratory analysis. On the other hand, this approach has the disadvantage that it presumes the SME knows everything there is to know about the data (in my experience this is rarely, if ever, the case), thereby precluding the possibility of serendipitously uncovering new relationships.
There is a third approach that is less reliant on SMEs and more empirically driven. This method assumes, as you have noted, that a relative handful of features are explanatorily important, i.e., they rank towards the top of relative variable importance. By focusing interactions on that subset of features, it becomes possible to derive a more manageably sized set of interactions, as long as there aren't too many important, candidate features.
So, how might one address your concern with pulling together an "ensemble" of metrics providing statistically useful information regarding the relationships between and among the features and their interactions. For me, this is entirely a function of the analyst's appetite, patience and time available for playing around with the results from an ensemble of differently specified models. A simple example of how this could work is OLS ANOVA with two features and one two-way interaction. Depending on the software, it is possible to specify that the sums of squares of the features and interaction be treated independently in entering the model -- some packages describe this as SS3 (for sums of squares type 3). One could also try SS1, SS2 or SS4 which invoke different orderings of the independent variables. By playing around with the permutations of these orderings, the "humanly interpretable" information can be deduced...that is, up to the limits of human cognition in juggling a large number of differing metrics. Ulrike Groemping has done the most work in reviewing and explaining the different statistical approaches to this kind of model evaluation of relative variable importance. She also has an R module which does a computationally intensive version of this called RELAIMPO. Check her papers out.
My personal solution is different from hers. Imagine you have summary performance statistics across a set of different model specifications returning metrics such as score, r-square, AIC, different specifications of SS types, and so on. Making sure they are all in a consistent unit of analysis would be a necessary step. Once that's done and with enough observed summary values, a PCA of these metrics becomes feasible. By scoring the input data based on the first PCA, then a summary ranking of the units under analysis is possible. I'm not sure I'm explaining this well enough but, having gotten it to work, it can be a useful tool for information reduction particularly in the case when one has a set of competing summary metrics. Pushback from other CV participants is expected here.
