What could cause ML predictions to be multimodal when the inputs are unimodal? 
This appears to be a systematic issue that occurs with the whole range of data points:

In general, what are some possible causes for discrepancies like this? I know it's impossible to pinpoint the exact cause without seeing my script and data, but I'm clueless as to how to even start debugging this.
I'm using XGBoost.
 A: This can happen for at least two different reasons.


*

*Your predictor has different distributions in the training and the prediction step. If your predictor values during training cluster in one part of space, you have many training response values corresponding to this part of predictor space. And of course, if your predictor values cluster in a different part of space during prediction, your predicted response will cluster around values in that part of predictor space.
So check whether the distributions of your predictors are similar.

*Now, let's assume the distributions of your predictors are similar. One reason for the effect you observe could be that the actual mean response to your predictors is multimodal, but the actual responses you see may be swamped by noise. If you have enough data, your model may pick up on the actual underlying multimodal signal, and your (point!) prediction may reflect this.
An example. Here is the actual mean response to a single predictor:

Let's assume we have 5,000 observations of the predictor, which are uniformly distributed in [-1,1]. We also have the corresponding values of the response. Because of noise, the response is unimodal (note that this histogram does not show the predictor any more!):

Since we have a lot of data, our model picks up on the multimodal relationship between the predictor and the response. Here is a plot of the predicted response against the predictor:

(Note how this is much wigglier than the actual relationship above. There is a lot of overfitting going on here. But that is not the main issue.)
If we now again take uniformly distributed predictor values and plot a histogram of the associated predicted responses, we get multimodality:

Bottom line: don't just look at histograms of responses. Plot observed and predicted responses against predictors. Also, consider looking not only at point predictions, but at full predicted densities, e.g., by sampling from them and plotting histograms from these resamples. If you do this for predictor values that are similarly distributed as in your training data, the resulting histograms should look similar to histograms of your training responses.
R code:
mean_response <- function(predictor) -.5*predictor^4+.4*predictor^2+.1
predictor_plot <- seq(-1,1,by=.01)

plot(predictor_plot,mean_response(predictor_plot),type="l",xlab="Predictor",ylab="Mean response")

set.seed(1)
predictor_train <- runif(5e3,-1,1)
response_train <- rnorm(length(predictor_train),mean_response(predictor_train),.1)
hist(response_train)

library(randomForest)
model <- randomForest(response_train~predictor_train)
predictions <- predict(model,newdata=data.frame(predictor_train=predictor_plot))
plot(predictor_plot,predictions,type="l")
hist(predictions)

