Investigating differences between populations Say we have a sample from two populations: A and B. Let's assume these populations are made of individuals and we choose to describe individuals in terms of features. Some of these features are categorical (e.g. do they drive to work?) and some are numerical (e.g. their height). Let's call these features: $X_1 \ldots X_n$.  We collect hundreds of these features (e.g. n = 200), let's assume for simplicity, with no error nor noise across all individuals.
We hypothesize the two populations are different. Our goal is to answer the following two questions:


*

*Are they actually significantly different?

*What is significantly different between them?


Methods like decision trees (e.g. random forests) and linear regression analysis can help. For example, one could look at feature importance in random forests or the fitted coefficients in linear regression to understand what may distinguish these groups, and explore relationships between features and populations. 
Before I go down this route, I want to get a sense of my options here, what's good and modern vs bad practice. Please note that my goal isn't prediction per se, but testing and finding any significant differences between the groups.
What are some principled approaches to address this problem? 
Here are some concerns I have:


*

*Methods like linear regression analysis may not fully answer (2), right? E.g. a single fit can help find some differences, but not all significant differences. For example, multi-collinearity may prevent us from finding how all features vary across groups (at least on a single fit). For the same reason, I would expect ANOVA can't provide a full answer to (2) either.

*It's not entirely clear how a predictive approach would answer (1). For example, what classification/prediction loss function should we minimize? And how do we test whether or not the groups are significantly different once we have a fit? Finally, I worry that the answer I get to (1) may depend on the particular set of classification models I use.
 A: Let's think the problem as follows.
Say $X=(X_1,X_2,..X_n)$ and $Y$ is a binary variable standing for the population : $Y=0$ means first population, $Y=1$ means second population. The null hypothesis can be expressed in several equivalent ways:


*

*$H_0$: the populations are the same

*$H_0$:
the distribution of $X$ given $Y=0$ is the same as the distribution of $X$ given $Y=1$

*$H_0$: $X$ and $Y$ are independent

*$H_0$: for any function $f$ into $\{0,1\}$, $f(X)$ and $Y$ are independent


I don't know much about random forests, but they may be thought as an all purpose predictor that avoids over-fitting. If we idealize them quite a bit: it is something capable of detecting any kind of relationship between $Y$ and any kind of features $X$ without over-fitting.
It is possible to try something based on this. Split the original dataset into a training set and a test set. Then:


*

*train a random forest $f$ that predicts $Y$ from $X$ on the training set.

*make a simple chi-squared independence test (with risk $\alpha$) between $f(X)$ and $Y$ on the test set


This test is quite conservative. If the random forest is a poor method, at worst outputting a dumb $f(X)$, then it will reject $H_0$ with a probability less than $\alpha$ anyway (when $H_0$ is true). The over-fitting would not even be a problem since we use a test and a training set. However, the power of the test directly depends on the intelligence of the random forest method (or any predictor used).
Note that you can use several possible predictors: like plain old logistic regression first, then logistic regression with some cross features, then a few decision trees, then a random forest... But if you do so you should adjust $\alpha$ to the number of tests to avoid "false discoveries". See: Alpha adjustment for multiple testing
A: You don't say how many features are available in the data. Few, many, massive? Can we assume they are the same features between populations, all measured using the same tools, methods and modalities? If not, then you have a bigger problem where an errors-in-variables measurement model might work.
@benoitsanchez appears to have answered question #1). 
Wrt #2), I'm not sure RFs can help. By using a more formal model such as one-way ANOVA applied to one feature at a time, a test of the difference between populations for features can be developed. By summarizing the results of those tests, based on the magnitude of the test as well as its significance, a descriptive profile of how the populations differ across features becomes possible. This is an admittedly ad hoc and heuristic solution which may not be rigorous enough for your tastes, preferences and training. 
Not being good at Latex-type notation, let me simply describe how these tests might work: first, construct some kind of macro loop which passes all features through, one feature at a time. With each pass of the loop, the new feature becomes the target or DV with X consisting of a dummy variable for population as well as any control variables that are appropriate. Make sure that the same controls are used for each feature as well as that the underlying data is exactly the same for all ANOVAs, eliminating variation attributable to the vicissitudes of finite data samples. Aggregate the F-test values for the dummy variable for each feature. This will provide a standardized metric enabling comparison across features. F-tests are preferable to fitted betas since betas are not standardized, being expressed in the unit and std devs of each individual feature. 
Your last comment, "I worry that the answer I get to (1) may depend on the particular set of classification/regression models that I use," is always true. The answers are quite likely to vary as a function of the model(s) used. It is also an expression of a commonly observed malaise among the more strongly theoretical and classically trained statisticians who aren't comfortable with or have trouble acknowledging the non-deterministic nature of applied statistical modeling. An excellent antidote for these symptoms is Efron and Hastie's recent book Computer Age Statistical Inference. They bring statistical modeling into the 21st c, an age of data science and machine learning, by candidly acknowledging the iterative, approximating, heuristic nature of all models possessing an error term. One doesn't have to be a Bayesian to recognize the truth inherent in this observation. Their's is a refreshing perspective that differs from the rigid determinism of classical, 20th c statistical practice which threw up its hands when, e.g., a cross-products matrix wouldn't invert and/or some pedantic model assumption wasn't met. 
