  Motivated by models for multiway comparison data, we consider the problem of
estimating a coordinate-wise isotonic function on the domain $[0, 1]^d$ from
noisy observations collected on a uniform lattice, but where the design points
have been permuted along each dimension. While the univariate and bivariate
versions of this problem have received significant attention, our focus is on
the multivariate case $d \geq 3$. We study both the minimax risk of estimation
(in empirical $L_2$ loss) and the fundamental limits of adaptation (quantified
by the adaptivity index) to a family of piecewise constant functions. We
provide a computationally efficient Mirsky partition estimator that is minimax
optimal while also achieving the smallest adaptivity index possible for
polynomial time procedures. Thus, from a worst-case perspective and in sharp
contrast to the bivariate case, the latent permutations in the model do not
introduce significant computational difficulties over and above vanilla
isotonic regression. On the other hand, the fundamental limits of adaptation
are significantly different with and without unknown permutations: Assuming a
hardness conjecture from average-case complexity theory, a
statistical-computational gap manifests in the former case. In a complementary
direction, we show that natural modifications of existing estimators fail to
satisfy at least one of the desiderata of optimal worst-case statistical
performance, computational efficiency, and fast adaptation. Along the way to
showing our results, we improve adaptation results in the special case $d = 2$
and establish some properties of estimators for vanilla isotonic regression,
both of which may be of independent interest.
