  This paper studies the Gaussian and bootstrap approximations for the
probabilities of a non-degenerate U-statistic belonging to the hyperrectangles
in $\mathbb{R}^d$ when the dimension $d$ is large. A two-step Gaussian
approximation procedure that does not impose structural assumptions on the data
distribution is proposed. Subject to mild moment conditions on the kernel, we
establish the explicit rate of convergence uniformly in the class of all
hyperrectangles in $\mathbb{R}^d$ that decays polynomially in sample size for a
high-dimensional scaling limit, where the dimension can be much larger than the
sample size. We also provide computable approximation methods for the quantiles
of the maxima of centered U-statistics. Specifically, we provide a unified
perspective for the empirical bootstrap, the randomly reweighted bootstrap, and
the Gaussian multiplier bootstrap with the jackknife estimator of covariance
matrix as randomly reweighted quadratic forms and we establish their validity.
We show that all three methods are inferentially first-order equivalent for
high-dimensional U-statistics in the sense that they achieve the same uniform
rate of convergence over all $d$-dimensional hyperrectangles. In particular,
they are asymptotically valid when the dimension $d$ can be as large as
$O(e^{n^c})$ for some constant $c \in (0,1/7)$.
  (Full abstract can be found in the paper.)
