A couple of thoughts on this problem which might be of interest: your predictor can be interpreted as the Bayes classifier for a Gaussian mixture model. For example, if you take $r \in \mathbf{R}^d, Q \succ 0$, then
\begin{align}
y_j(x) &= \frac{\omega_j \mathcal{N} \left( x | \mu_j, Q^{-1} \right)}{\sum_{k=1}^n \omega_k \mathcal{N} \left( x | \mu_k, Q^{-1} \right)} \\
\text{where} \quad\omega_k &= \exp \left( b_k + \frac{1}{2} |a_k + r|_{Q^{-1}}^2 \right) \\
\mu_k &= Q^{-1} \left( a_k + r \right).
\end{align}
Your quantity of interest can be expressed as "given a uniform sample from the hypercube, what is the probability that it came from the $j^{\text{th}}$ component of this GMM?". This naturally suggests a simple importance sampling scheme which puts most of its mass near the $j^{\text{th}}$ mode, e.g.
\begin{align}
q_j (x) = \alpha_0 \cdot \mathbb{I} \left[ x \in [0,1]^d \right] + \alpha_1 \cdot \mathcal{N} \left( x | \mu_j, \gamma \cdot Q^{-1} \right)
\end{align}
where $(\alpha_0, \alpha_1, \gamma)$ are tunable parameters. One could also consider importance distributions based on the full $n$-component mixture, e.g.
\begin{align}
q^n (x) = \alpha_0 \cdot \mathbb{I} \left[ x \in [0,1]^d \right] + \sum_{k = 1}^n \alpha_k \mathcal{N} \left( x | \mu_k, \gamma \cdot Q^{-1} \right),
\end{align}
which induces additional tuning parameters, which has its costs and benefits. Of course, there is the complication of setting $(r, Q)$ as well. My guess is that if one can set $(r, Q)$ such that the support of the implied GMM overlaps well with the hypercube as a whole (by some measure), this will give well-behaved weights for a derived importance sampling scheme.
My suspicion here is that if one can find a good way of setting $(r, Q)$, then life is not too hard: one can use the above construction to parametrise a family of importance densities, and then apply standard adaptive importance sampling techniques, possibly also using some extra variance reduction tricks like QMC / stratification and control variates to sweeten the deal.
I hope this is at least somewhat useful. I am a little frustrated that I lack intuition for how to set $(r, Q)$ automatically, as once that is done, the problem is relatively friendly. One idea would be to try to tune $(r, Q)$ such that
\begin{align}
\sum_{k = 1}^K w_k \textbf{Var}_{\mathcal{N} \left( x | \mu_k, Q^{-1} \right)} \left( y_j (x) \cdot \mathbb{I} \left[ x \in [0,1]^d \right] \right)
\end{align}
is minimal, though it should be emphasised that this is just a heuristic.
