summaries: This paper studies a problem called lattice partition recovery, where piecewise constant signals over a d-dimensional lattice are corrupted by noise and require to be recovered. A method called dyadic classification and regression tree (DCART) is studied. It is shown that DCART achieves one-sided consistency, but suffers from over-partitioning. To achieve partition recovery, this paper proposes an additional merging stage, where dyadic rectangles that are close are merged. This is shown to achieve partition recovery under mild conditions. A minimax lower bound is also given, and numerical experiments are shown to support the methodology. This paper considers a piecewise constant signal on a $d$-dimensional lattice. The observed signal is corrupted by additive Gaussian noise. The task is to recover the partition of the lattice such that the signal within a partition is constant. The paper proposes an algorithm that builds upon the scalable DCART algorithm.  The contributions of the paper are threefold:  (i)	The authors first analyze the performance of the vanilla DCART algorithm to recover the partitions. They prove the one-sided consistency of the vanilla DCART algorithm. One-sided consistency guarantees that the recovered partitions have large subsets which have constant signal, but a partition of the ground truth signal could potentially be over-partitioned by the DCART algorithm. They demonstrate by an example that it is indeed the case.  (ii)	In order to address the issue of over-partitioning, the authors propose a two-step algorithm where the first step is the vanilla DCART algorithm while restricting to partitions with size above a threshold. In the second step, the partitions which are close and whose signals have a small difference are merged. They prove that this estimator is consistent.  (iii)	They prove that the proposed 2 step estimator is minimax optimal if the number of partitions is constant. Finally, they demonstrate the efficacy of their algorithm by simulations.  The authors consider the problem of recovering the partition underlying a piecewise constant function over a $d$-dimensional lattice, _i.e._, a signal with support $L_{d,n} = \lbrace 1, \dots, n\rbrace^d$. That is, we wish to recover the structure such that within each rectangular partition of the lattice, the function $\theta^* \in \mathbb{R}^{L_{d,n}}$ is piecewise constant, from noisy measurements ${y_i = \theta_{i}^* + \epsilon_i}$ for $i \in L_{d,n}$ where $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$ i.i.d. A straightforward formulation would be NP-hard, so they the consider the dyadic CART (DCART) algorithm, which simplifies the search space by restricting it to partitions that can be formed by repeated bisections of partitions, starting with a single partition containing all $N = n^d$ points in the lattice.  The authors provide a number of bounds of the form $\sigma^2 k_{\mathrm{dyad}}(\theta^*) \log(N) / \kappa^2$, where $\kappa$ is the minimum jump size (difference between signal in adjacent partitions) of $\theta^*$. First, they show that using DCART, the number of lattice points in recovered partitions where the true signal is not actually piecewise constant scales at most with the above expression ("one-sided consistency"). They then propose a two-step estimator that adds a step merging rectangles in recovered partitions whose measurements are "close" to each other, and show that under certain assumptions, the Hausdorff distance between the recovered (estimated) partition and the actual partition (accounting for the need to merge rectangles to assure definitional uniqueness) scales at most with the above expression. Finally, they provide a simple example to demonstrate an order-matching lower bound. In this work the authors consider the problem of learning the partition of a piecewise constant signal supported on a $d$-dimensional lattice $\{1, \ldots, n\}^d$.  They assume access to a noisy signal $y = \theta^\ast + \epsilon$, where $\epsilon \sim \mathcal{N}(0, \sigma^2 I)$ is a Gaussian noise vector and $\theta^\ast$ is the true signal that is assumed to be piecewise constant over a rectangular partition of the latice.  The goal is to recover the underlying partition.  This work focuses on the performance of (variants of) the Dyadic CART (DCART) algorithm for the task of partition recovery.  Prior work [1] has shown that the DCART estimator achieves the optimal error for the task of recovering the true signal $\theta^\ast$.  The authors argue that the partition produced by the vanilla DCART algorithm may contain a large number of rectangles.  However, they show that inside each rectangle of the partition of the vanilla DCART, the signal is indeed almost constant.  The second result of this work proposes simple modifications for DCART so that the corresponding partition is close to the true underlying partition.  In particular they enforce a lower bound on the size of the rectangles of DCART and, at the end, merge cells that are close and have similar value.  They show that the modified DCART algorithm outputs a partition at distance roughly $\sigma^2 k^\ast \log(N) / \kappa^2$, where $N = n^d$ is the size of the grid, $\kappa$ is the minimum signal difference between contiguous cells of the true signal, and $k^\ast$ is the minimum number of dyadic cells required to create a partition consistent with the true signal.  The authors show a lower bound of $\sigma^2 \log(N)/\kappa^2$ for the error. Therefore, their upper and lower bounds differ by a factor of $k^\ast$.  Finally, the authors are able to remove the extraneous factor $k^\ast$ under the additional assumption that the underlying partition is not very dense, i.e., around every rectangle of the true underlying partition there are not many other rectangles.   [1]: Sabyasachi Chatterjee and Subhajit Goswami Adaptive estimation of multivariate piecewise polynomials and bounded variation functions by optimal decision trees This article investigates the problem of the recovery of constancy regions of a noisy high dimensional piecewise constant signal. The authors proposed to use the DCART estimator to build up a rectangular partition and proved a weak consistency result. Then, they proposed a two steps estimator based on an appropriate merging of the rectangles that DCART output. For this variant (DCARTAM), they proved consistency: the estimator converges to the ground truth partition in the sense of Hausdorff (after an appropriate normalization by the number of rectangles in the partitions).  Finally, DCART and DCARTAM were compared to the TV-based estimator empirically using two metrics. The first one is the Hausdorff distance between the estimator and the ground truth partition, and the second one is the difference between the cardinality of the estimator and the cardinality of the ground truth partition. Visually, DCARTAM seems to be superior to DCART. Moreover, DCARTAM seems to be superior to the TV-based estimator using the second distance, yet it does not show a systematic superiority if the first distance was used.
abstractText: We study piece-wise constant signals corrupted by additive Gaussian noise over a d-dimensional lattice. Data of this form naturally arise in a host of applications, and the tasks of signal detection or testing, de-noising and estimation have been studied extensively in the statistical and signal processing literature. In this paper we consider instead the problem of partition recovery, i.e. of estimating the partition of the lattice induced by the constancy regions of the unknown signal, using the computationally-efficient dyadic classification and regression tree (DCART) methodology proposed by [14]. We prove that, under appropriate regularity conditions on the shape of the partition elements, a DCART-based procedure consistently estimates the underlying partition at a rate of order σ2k∗ log(N)/κ, where k∗ is the minimal number of rectangular sub-graphs obtained using recursive dyadic partitions supporting the signal partition, σ is the noise variance, κ is the minimal magnitude of the signal difference among contiguous elements of the partition and N is the size of the lattice. Furthermore, under stronger assumptions, our method attains a sharper estimation error of order σ log(N)/κ, independent of k∗, which we show to be minimax rate optimal. Our theoretical guarantees further extend to the partition estimator based on the optimal regression tree estimator (ORT) of [12] and to the one obtained through an NP-hard exhaustive search method. We corroborate our theoretical findings and the effectiveness of DCART for partition recovery in simulations.
authors.name: Oscar Hernan Madrid Padilla
authors.name: Yi Yu
references.authors: Louigi Addario-Berry Nicolas Broutin Luc Devroye Gábor Lugosi
references.title: On combinatorial testing problems
references.venue: The Annals of Statistics,
references.authors: Ery Arias-Castro Emmanuel J Candes Arnaud Durand
references.title: Detection of an anomalous cluster in a network
references.venue: The Annals of Statistics,
references.authors: Ery Arias-Castro Emmanuel J Candes Hannes Helgason Ofer Zeitouni
references.title: Searching for a trail of evidence in a maze
references.venue: The Annals of Statistics,
references.authors: Ery Arias-Castro Emmanuel J. Candès Yaniv Plan
references.title: Global testing under sparse alternatives: Anova, multiple comparisons and the higher criticism
references.venue: The Annals of Statistics,
references.authors: Ery Arias-Castro Yuchao Liu
references.title: Distribution-free detection of a submatrix
references.venue: Journal of Multivariate Analysis,
references.authors: JiaWang Bian Wen-Yan Lin Yasuyuki Matsushita Sai-Kit Yeung Tan-Dat Nguyen Ming- Ming Cheng
references.title: Gms: Grid-based motion statistics for fast, ultra-robust feature correspondence
references.venue: In Proceedings of the IEEE conference on computer vision and pattern recognition,
references.authors: Leo Breiman Jerome Friedman Charles J Stone Richard A Olshen
references.title: Classification and regression trees
references.venue: CRC press,
references.authors: Victor-Emmanuel Brunel
references.title: Adaptive estimation of convex polytopes and convex sets from noisy data
references.venue: Electronic Journal of Statistics,
references.authors: C. Butucea Y.I. Ingster I.A. Suslina
references.title: Sharp variable selection of a sparse submatrix in a high-dimensional noisy matrix
references.venue: ESAIM: Probability and Statistics,
references.authors: Cristina Butucea Yuri I. Ingster
references.title: Detection of a sparse submatrix of a high-dimensional noisy matrix
references.venue: Bernoulli, 19(5B):2652–2688,
references.authors: T. Tony Cai Tengyuan Liang Alexander Rakhlin
references.title: Computational and statistical boundaries for submatrix localization in a large noisy matrix
references.venue: The Annals of Statistics,
references.authors: Sabyasachi Chatterjee Subhajit Goswami
references.title: Adaptive estimation of multivariate piecewise polynomials and bounded variation functions by optimal decision trees
references.authors: Yudong Chen Jiaming Xu
references.title: Statistical-computational tradeoffs in planted problems and submatrix localization with a growing number of clusters and submatrices
references.venue: Journal of Machine Learning Research,
references.authors: David L Donoho
references.title: Cart and best-ortho-basis: a connection
references.venue: The Annals of statistics,
references.authors: Zhou Fan Leying Guan
references.title: Approximate `0-penalized estimation of piecewise-constant signals on graphs
references.venue: The Annals of Statistics,
references.authors: Evelina Fedorenko John Duncan Nancy Kanwisher
references.title: Broad domain generality in focal regions of frontal and parietal cortex
references.venue: Proceedings of the National Academy of Sciences,
references.authors: Chao Gao Yu Lu Zongming Ma Harrison H. Zhou
references.title: Optimal estimation and completion of matrices with biclustering structures
references.venue: Journal of Machine Learning Research,
references.authors: Bruce Hajek Yihong Wu Jiaming Xu
references.title: Submatrix localization via message passing
references.venue: Journal of Machine Learning Research,
references.authors: Qiyang Han
references.title: Global empirical risk minimizers with "shape constraints" are rate optimal in general dimensions
references.venue: arXiv preprint arXiv:1905.12823,
references.authors: Abhishek Kaul
references.title: Segmentation of high dimensional means over multi-dimensional change points and connections to regression trees
references.venue: arXiv preprint arXiv:2105.10017,
references.authors: Mladen Kolar Sivaraman Balakrishnan Alessandro Rinaldo Aarti Singh
references.title: Minimax localization of structural information in large noisy matrices
references.venue: Advances in Neural Information Processing Systems
references.authors: Aleksandr P Korostelev Aleksandr B Cybakov
references.title: Asymptotically minimax image reconstruction problems
references.venue: Sonderforschungsbereich 123,
references.authors: Andrew Lang Aaron Carass Peter A Calabresi Howard S Ying Jerry L Prince
references.title: An adaptive grid for graph-based segmentation in retinal oct
references.venue: In Medical Imaging 2014: Image Processing,
references.authors: Yuchao Liu Ery Arias-Castro
references.title: A multiscale scan statistic for adaptive submatrix localization
references.venue: arXiv preprint arXiv:1906.08884,
references.authors: Zongming Ma Yihong Wu
references.title: Computational barriers in minimax submatrix detection
references.venue: The Annals of Statistics,
references.authors: Vincent Roullier Olivier Lézoray Vinh-Thong Ta Abderrahim Elmoataz
references.title: Multi-resolution graph-based analysis of histopathological whole slide images: Application to mitotic cell extraction and visualization
references.venue: Computerized Medical Imaging and Graphics,
references.authors: Leonid I Rudin Stanley Osher Emad Fatemi
references.title: Nonlinear total variation based noise removal algorithms
references.venue: Physica D: nonlinear phenomena,
references.authors: Andrey A. Shabalin Victor J. Weigman Charles M. Perou Andrew B. Nobel
references.title: Finding large average submatrices in high dimensional data
references.venue: The Annals of Applied Statistics,
references.authors: James Sharpnack Aarti Singh Alessandro Rinaldo
references.title: Sparsistency of the edge lasso over graphs
references.venue: In Artificial Intelligence and Statistics,
references.authors: Jonathan R Stroud Michael L Stein Shaun Lysen
references.title: Bayesian and maximum likelihood estimation for gaussian processes on an incomplete lattice
references.venue: Journal of computational and Graphical Statistics,
references.authors: Xing Sun Andrew B. Nobel
references.title: On the maximal size of large-average and ANOVA-fit submatrices in a Gaussian random matrix
references.authors: Wesley Tansey Oluwasanmi Koyejo Russell A Poldrack James G Scott
references.title: False discovery rate smoothing
references.venue: Journal of the American Statistical Association,
references.authors: Wesley Tansey James G Scott
references.title: A fast and flexible algorithm for the graph-fused lasso
references.venue: arXiv preprint arXiv:1505.06475,
references.authors: Ryan J Tibshirani Jonathan Taylor
references.title: The solution path of the generalized lasso
references.venue: The annals of statistics,
references.authors: Nicolas Verzelen Magalie Fromont Matthieu Lerasle Patricia Reynaud-Bouret
references.title: Optimal change-point detection and localization
references.venue: arXiv preprint arXiv:2010.11470,
references.authors: Daren Wang Yi Yu Alessandro Rinaldo
references.title: Univariate mean change point detection: Penalization, cusum and optimality
references.venue: Electronic Journal of Statistics,
references.authors: Timothy G Whiteside Andrew J Esparon Renée E Bartolo
references.title: A semi-automated approach for quantitative mapping of woody cover from historical time series aerial photography and satellite imagery
references.venue: Ecological Informatics,
references.authors: Sascha Wirges Tom Fischer Christoph Stiller Jesus Balado Frias
references.title: Object detection and classification in occupancy grid maps using deep convolutional networks
references.venue: In 2018 21st International Conference on Intelligent Transportation Systems (ITSC),
sections.heading: 1 Introduction
sections.text: Suppose we observe a noisy realization of a structured, piece-wise constant signal supported over a d-dimensional square lattice (or grid graph) Ld,n = {1, . . . , n}d. Data that can be modeled in this manner arise in several application areas, including in satellite imagery [e.g. 30, 37], computer vision [e.g. 6, 38], medical imaging [e.g. 26, 23], and neuroscience [e.g. 16, 32]. Our goal is to estimate the constancy regions of the underlying signal. Specifically, we assume that the data y ∈ RLd,n are such that, for each coordinate i ∈ Ld,n,
  yi = θ ∗ i + i, (1)
  where ( i, i ∈ Ld,n) are i.i.d. N (0, σ2) noise variables and the unknown signal θ∗ ∈ RLd,n is assumed to be piece-wise constant over an unknown rectangular partition of Ld,n. We define a subset R ⊂ Ld,n to be a rectangle if R = ∏d i=1[ai, bi], where [a, b] = {j ∈ Z : a ≤ j ≤ b}, a, b ∈ Z. A rectangular partition of Ld,n, P , is a collection of disjoint rectangles {Rl} ⊂ Ld,n, satisfying ∪R∈PR = Ld,n. To each vector in Ld,n, there corresponds a (possibly trivial) rectangular partition.
  35th Conference on Neural Information Processing Systems (NeurIPS 2021).
  Definition 1 A rectangular partition associated with a vector θ ∈ RLd,n is a rectangular partition {Rl}l∈[1,k] of Ld,n, such that θ takes on constant values over each Rl. For a vector θ ∈ RLd,n , we let k(θ) be the smallest positive integer such that there exists a rectangular partition with k(θ) elements and associated with θ.
  In this paper, we are interested in recovering a rectangular partition associated with the signal θ∗ in (1). A complication immediately arises when d ≥ 2: the rectangular partition associated with a given θ is not necessarily unique. This fact is illustrated in Figure 1, where the left-most plot depicts the lattice supported vector θ, consisting of a rectangle of elevated value (in grey) against a background (in white). For such θ, we show three possible rectangular partitions, each of which comprised of five rectangles (the second, third and fourth plots). In fact, the partition recovery problem is well defined, as long as we consider coarser partitions comprised by unions of adjacent rectangles instead of individual rectangles: see Definition 2 below for details. We remark that this issue does not occur in the univariate (d = 1) case, for which the partition recovery task has been thoroughly studied in the change point literature; see section 1.3 below. Thus, we assume that d ≥ 2. For the purpose of estimating the rectangular partition associated with θ∗ (or, more precisely, its unique coarsening as formalized in Definition 2), we resort to the dyadic classification and regression tree (DCART) algorithm of [14]. This is a polynomial-time decision-tree-based algorithm developed for de-noising purposes for signals over lattices, and is a variant of the classification and regression trees (CART) [7]. See Section 1.1 below for a description of DCART. The optimal regression trees (ORT) estimator, recently proposed in [12], further builds upon DCART and delivers sharp theoretical guarantees for signal estimation while retaining good computational properties – though we should mention that in our experiments we have found DCART to be significantly faster. Both DCART and its more sophisticated version ORT can be seen as approximations to the NP-hard estimator
  θ1 = arg minθ∈RLd,n { 2−1‖y − θ‖2 + λk(θ) } , (2)
  where k(θ) is given in Definition 1, ‖ · ‖ is the vector (or Euclidean) `2-norm and λ > 0 a tuning parameter. DCART modifies the above, impractical optimization problem by restricting only to dyadic rectangular partitions. This leads to significant gains in computational efficiency without sacrificing on the statistical performance. Indeed, decision-tree-based algorithms have been shown to be optimal under various settings for the purpose of signal estimation; see [12, 20]. In this paper, we further demonstrate their effectiveness for the different task of partition recovery. In particular, we show how simple modifications of the DCART (or ORT) estimator yield practicable procedures for partition recovery with good theoretical guarantees and derive novel localization rates.
  Note that, there is a wide array of applications focusing on detecting the regions rather than estimating the background signals, especially in surveillance and environment monitoring. Our work is motivated by all the applications/problems considered in the large literature on biclustering, where the underlying signal is assumed to be piecewise constant. Estimating the boundary of the partition is the most refined and difficult task in these settings. Thus, any of the many scenarios in which biclustering is relevant can be used to motive our task. An analogous observation holds also for the more general problem of identifying an anomalous cluster (sub-graph) in a network, a problem that has been tackled (for testing purposes only) by [4], the reference therein provide numerous examples of applications. On a high-level, the relationship between the partition and signal recoveries can be thought of the relationship between the estimation consistency and support consistency in a high-dimensional linear regression problems. They can be done by almost identical algorithms but the theoretical results rely on different sets of conditions.
  The paper is organized as follows. In the rest of this section we formalize the problem settings and the task of partition recovery, and describe the DCART procedure. We further summarize our main findings and discuss related literature. Section 2 contains our main results about one- and
  two-sided consistency of DCART and its modification. In Section 2.3 we derive a minimax lower bound stemming from the case of one rectangular region of elevated signal against pure background. Illustrative simulations corroborating our findings can be found in Section 3. The Supplementary Material contains the proofs.
  Notation We set N = nd, the size of the lattice Ld,n, where we recall that d ≥ 2 is assumed fixed throughout. For any integer m ∈ N∗, let [m] = [1,m]. Given a rectangular partition Π of Ld,n, let S(Π) be the linear subspace of RLd,n consisting of vectors with constant values on each rectangle in Π and let OS(Π)(·) be the orthogonal projection onto S(Π). For any R ⊂ Ld,n and θ ∈ RLd,n , let θ̄R = |R|−1 ∑ i∈R θi, where | · | is the cardinality of a set. Two rectangles R1, R2 ∈ Π are said to be adjacent if there exists l ∈ [d] such that R1 and R2 share a boundary along el and one is a subset of the other in the hyperplane defined by el, the lth standard basis vector in Rd. See Definition S1 for a rigorous definition. This concept of adjacency is specifically tailored to – and in fact only valid for – dyadic (and hierarchical, in the sense specified by [12]) rectangular partitions, which are most relevant for this paper. For any subsets A,B ⊂ Ld,n, define dist(A,B) = mina∈A,b∈B ‖a − b‖. Throughout this paper, we will use the `2-norm as the vector norm.
sections.heading: 1.1 Problem setup
sections.text: We begin by introducing two key parameters for the model specified in (1) and a well-defined notion of rectangular partition induced by θ∗.
  Definition 2 (Model parameters, induced partitions) Let θ∗ as in (1) and {R∗j}j∈[m] be a rectangular partition of Ld,n associated with θ∗. Consider the graphG∗ = (E∗, V ∗), where V ∗ = [m] and E∗ = {(i, j) : θ̄∗R∗i = θ̄ ∗ R∗j , R∗i and R ∗ j are adjacent}. Let {C∗l }l∈[L] be all connected components of G∗ and define Λ∗ = Λ∗(θ∗) = {∪j∈C∗1R ∗ j , . . . ,∪j∈C∗LR ∗ j} as the partition (not necessarily rectangular) induced by θ∗. We say that the union of rectangles ∪j∈C∗sR ∗ j and ∪j∈C∗t R ∗ j , s, t ∈ [L], s 6= t, are adjacent, if and only if there exists (i, j) ∈ C∗s × C∗t such that R∗i and R∗j are adjacent. Let κ and ∆ be the minimum jump size and minimal rectangle size, respectively, formally defined as
  κ = min a∈A,b∈B, A,B∈Λ∗, θ∗a 6=θ∗b |θ∗a − θ∗b | and ∆ = min j∈[m] |R∗j |.
  It is important to emphasize the difference between a partition associated with θ∗, as described in Definition 1, which may not be unique, and the partition Λ∗ induced by θ∗ of Definition 2, which is instead unique and thus describes a well-defined functional of θ∗. The parameters κ and ∆ capture two complementary aspects of the intrinsic difficulty of the problem of estimating Λ∗; intuitively, one would expect the partition recovery task to be more difficult when κ and ∆ are small (and σ is large). Below, we will prove rigorously that this intuition is indeed correct. When d = 1, both parameters, along with σ, have in fact been shown to fully characterize the change point localization task: see, e.g., [36, 35].
  The partition recovery task can therefore be formulated as that of constructing an estimator Λ̂ of Λ∗, the induced partition of θ∗, such that, as the sample size N grows unbounded and with probability tending to one,
  |Λ̂| = |Λ∗| and ∆−1dHaus(Λ̂,Λ∗) = ∆−1 max A∈Λ∗ min B∈Λ̂ |A4B| → 0, (3)
  where A4B is the symmetric difference between A and B. We refer to dHaus(Λ̂,Λ∗) as the localization error for the partition recovery problem.
  The dyadic classification and regression trees (DCART) estimator. In order to produce a computationally efficient estimator of Λ∗ satisfying the consistency requirements (3), we deploy the DCART procedure [14], which can be viewed as an approximate solution to the problem in (2). Instead of optimizing over all vectors in RLd,n , DCART minimizes the objective function only over vectors associated with a dyadic rectangular partition, which is defined as follows. LetR = ∏ i∈[d][ai, bi] ⊂ Ld,n be a rectangle. A dyadic split of R chooses a coordinate j ∈ [d], l the middle point of [aj , bj ], and
  splits R into R1 = ∏
  i∈[j−1]
  [ai, bi]× [aj , l]× ∏
  i∈[j+1,d]
  [ai, bi] and R2 = ∏
  i∈[j−1]
  [ai, bi]× [l+1, bj ]× ∏
  i∈[j+1,d]
  [ai, bi],
  with [0] = [c2, c1] = ∅, c2 > c1. Assuming that n is a power of 2, starting from Ld,n itself, we proceed iteratively as follows. Given the partition {Ru}u∈[k], one chooses a rectangle Ru and performs a dyadic split on Ru that leads to the largest reduction in the objective function. Any partition constructed through a sequence of such steps is called a dyadic rectangular partition. With a pre-specified λ > 0, the DCART estimator is
  θ̃ = OS(Π̃)(y), where Π̃ ∈ arg min Π∈Pdyadic,d,n
  { 2−1‖y −OS(Π)(y)‖2 + λ|Π| } , (4)
  where Pdyadic,d,n is the set of all dyadic rectangular partitions of Ld,n. As shown in [14] and [12], the DCART estimator can be obtained via dynamic programming with a computational cost of O(N). Given any solution to (4), a natural (though, as we will see, sub-optimal) estimator of the induced partition of θ∗ is Π̃, the partition associated with the resulting DCART estimator θ̃. Importantly, by the property of DCART, and using the fact that the Gaussian errors have a Lebesgue density, Π̃ is in fact a dyadic-rectangular partition and is unique with probability one, and thus the resulting estimator is well-defined. (Equivalently, the partition associated with θ̃ and the one induced by θ̃ coincide.)
sections.heading: 1.2 Summary of our results
sections.text: We briefly summarize the contributions made in this paper.
  One-sided consistency of DCART. Though DCART is known to be a minimax rate-optimal estimator of θ∗ [12], for the task of partition recovery its associated partition Π̃ has sub-optimal performance. Indeed, due to the nature of the procedure, it is easy to construct cases in which the DCART over-partitions. See Figure 2. In these situations, DCART falls short with respect to the target conditions for consistency described in (3). Nonetheless, it is possible to prove a weaker one-sided consistency guarantee, in the sense that every resulting DCART rectangle is almost constant. In detail, let R = {Rl}l∈[m̃] be the rectangular partition defined by Π̃ in (4). Then, we show in Section 2.1 that, for any Ri ∈ R, there exists Si ⊂ Ri such that θ∗t = θ∗u, u, t ∈ Si, and ∑ i∈[m̃] |Ri\Si| . κ−2σ2kdyad(θ∗) log(N). Throughout, the quantity kdyad(θ∗) refers to the smallest positive integer k such that there is a k-dyadic-rectangular-partition of Ld,n associated with θ∗.
  Two-sided consistency of DCART: A two-step estimator. In order to resolve the unavoidable overpartitioning issue with the naive DCART partition estimator and in order to prevent the occurrence of spurious clusters, we develop a more sophisticated two-step procedure. In the first step we use a variant of DCART that discourages the creation of rectangles of small volumes. In the second step, we apply a pruning algorithm merging rectangles when their values are similar and the rectangles are not far apart. With probability tending to one as N → ∞, the final output Λ̂ satisfies (3) with dHaus(Λ̂,Λ
  ∗) ≤ κ−2σ2kdyad(θ∗) log(N). This result is the first of its kind in the setting of lattice with arbitrary dimension d ≥ 2. This is shown in Section 2.2. Optimality: A regular boundary case. In Section 2.3, we consider the special case in which, for each rectangle in the rectangular partitions induced by θ∗ only has O(1)-many rectangles within distance of order σ2κ−2kdyad(θ∗) log(N). While more restrictive than the scenarios we study in Sections 2.1 and 2.2, this setting is broader than the ones adopted in the cluster detection literature [e.g. 2, 1]. In this case, with probability approaching one as N →∞, the estimator Λ̂ satisfies (3) and dHaus(Λ̂,Λ
  ∗) ≤ κ−2σ2 log(N). This error rate is shown to be minimax optimal, with a supporting minimax lower bound result in Proposition 3.
sections.heading: 1.3 Related and relevant literature
sections.text: The problem at hand is closely related to several recent research threads involving detection and estimation of a structured signal. When d = 2, our settings can be viewed as a generalization of those used for the purpose of biclustering, i.e. detection and estimation of sub-matrices. Though
  relatively recent, the literature on this topic is extensive, and the problem has been largely solved, both theoretically and methodologically. See, e.g., [28], [21], [10], [25], [31], [24], [5], [11], [9], [17], [18], [13] and [28].
  In the more general settings postulating a structured signal supported over a graph (including the grid graph), sharp results for the detection problem of testing the existence of a sub-graph or cluster in which the signal is different from the background are available in the literature: see, [3], [2], [1]. Concerning the estimation problem, [34], [29], [12], [15] and others, focused on de-noising the data and upper-bounding ‖θ̂ − θ∗‖∗, where θ̂ is an estimator of θ∗ and ‖ · ‖∗ is some vector norm. In yet another stream of work [e.g. 19, 8, 22] concerned with empirical risk minimization, the problem is usually formulated as identifying a single subset. More discussions can be found in Section S1.
  What sets our contributions apart from those in the literature referenced above, which have primarily targeted detection and signal estimation, is the focus on the arguably different task of partition recovery. As a result, the estimation bounds we obtain are, to the best of our knowledge, novel as they do not stem directly from the existing results.
  It is also important to mention how the partition recovery task can be cast as a univariate change point localization problem. Indeed, when d = 1, the two coincide; see [36, 35]. However, the case of d ≥ 2 becomes significantly more challenging due to the lack of a total ordering over the lattice. Consequently, our results imply also novel localization rates for change point analysis in multivariate settings.
sections.heading: 2 Consistency rates for the partition recovery problem
sections.text: In this section, we investigate the theoretical properties of DCART and of a two-step estimator also based on DCART for partition recovery. We remark that instead of DCART, it is possible to deploy the ORT estimator [12] or the NP-hard estimator (2) in our algorithms. Our theoretical results still hold by simply replacing the term kdyad(θ∗), in both the upper bound and the choice of tuning parameters, with the smallest k such that there is a k-hierarchical-rectangular-partition (khier(θ∗)) or k-rectangular-partition (k(θ∗)) of Ld,n associated with θ∗, respectively. Thus, using these more complicated methodologies that scan over larger classes of rectangular partitions will result in smaller upper bounds in Theorems 1, 2 and 4. See [12] for details about the relationship of k(θ∗), kdyad(θ∗) and khier(θ∗).
sections.heading: 2.1 One-sided consistency: DCART
sections.text: As illustrated in Figure 2, the DCART procedure will produce too fine a partition in many situations, even if the signal is directly observed (i.e., there is no noise). Thus, the naive partition estimator based on the constancy regions of the DCART estimator θ̃ as in (4) will inevitably suffer from the same drawback. Nonetheless, it is still possible to demonstrate a one-sided type of accuracy and even consistency for such a simple and computationally-efficient estimator. Specifically, in our next result we show that in every dyadic rectangle supporting the DCART estimator, θ∗ has almost constant mean. The reverse does not hold however, as there is no guarantee that every rectangle in the partition induced by the true signal θ∗ is mostly covered by one dyadic DCART rectangle.
  Theorem 1 Suppose that the data satisfy (1) and that θ̃ is the DCART estimator (4) obtained with tuning parameter λ = Cσ2 log(N), where C > 0 is a sufficiently large absolute constant. Let {Rj}j∈[k(θ̃)] be the associated partition. For any j ∈ [k(θ̃)], let Sj ⊂ Rj be the largest subset of Rj such that θ∗ is constant on Sj . Then there exist absolute constants C1, C2, C3, C4, C5 > 0 such that, with probability at least 1− C1 exp{−C2 log(N)}, the following hold:
  • global one-sided consistency:∑ j∈[k(θ̃)] |Rj\Sj | ≤ C3κ−2σ2kdyad(θ∗) log(N); (5)
  • local one-sided consistency: for any j ∈ [k(θ̃)], if Rj \ Sj 6= ∅, then
  |Rj\Sj | ≤ C4κ−2j σ 2kdyad(θ ∗ Rj ) log(N), (6)
  where κj = mins,t∈Rj : θ∗s 6=θ∗t |θ ∗ s − θ∗t |; and
  • control on over-partitioning:
  k(θ̃) ≤ 2kdyad(θ∗) + C5. (7)
  We remark that k(θ̃) = kdyad(θ̃) due to the construction of θ̃. Theorem 1 consists of three results. We have mentioned the over-partitioning issue of DCART. The bound (7) shows that the over-partitioning is upper bounded, in the sense that the size of the partition induced by DCART is in fact of the same order of the size of the dyadic rectangular partition associated with θ∗.
  For each resulting rectangle Rj , (6) shows that it is almost constant, in the sense that if the signal possesses different values in Rj , then Rj includes a subset Sj which has constant signal value and the size |Rj \ Sj | is upper bounded by κ−2j σ2kdyad(θ∗Rj ) log(N), where κj is the smallest jump size within Rj . We note that since kdyad(θ∗Rj ) ≤ kdyad(θ
  ∗), if kdyad(θ∗) is assumed to be a constant as in the cluster detection literature [e.g. 4], then for general d ∈ N∗, (6) has the same estimation error rate as that in the change point detection literature [e.g. 36, 35].
  The result in (6) provides an individual recovery error, with the individual jump size κj , while paying the price kdyad(θ∗Rj ). In (5) we show that globally, when we add up the errors in all resulting rectangles, the overall recovery error is of order κ−2σ2kdyad(θ∗) log(N). When d = 1, [35] shows the minimax rate of the L1-Wasserstein distance between the vectors of true change points and of the change point estimators is of order κ−2σ2K, where K is the true number of change points. Comparing with this result, (5) can be seen as delivering a “one-sided” nearly optimal rate, saving for a logarithmic factor.
  In the change point localization case, i.e. when d = 1, one can show that k(θ∗Rj ) ≤ 3 [see, e.g. 36] just assuming mild conditions on ∆. However, as soon as d ≥ 2 this is no longer the case. As an illustration, consider the left plot of Figure 2, where the whole rectangle is Rj and the white one in the middle is Sj . Without further constraints on each component, having only conditions on ∆ will not prevent a very fragmented boundary, which can increase the term kdyad(θ∗Rj ) in (6).
sections.heading: 2.2 Two-sided consistency: A two-step estimator
sections.text: As we have seen in Section 2.1, despite having the penalty term λ|Π| to penalize over-partitioning in the objective function (4), since the optimization of DCART only restricts to all dyadic partitions, the naive DCART estimator still suffers from over-partitioning. To address this issue, we propose a two-step estimator which builds upon DCART and merges the corresponding rectangles if they are close enough and if their estimated means are similar. The procedure will not only be guaranteed to return, with high probability, the correct number of rectangles in the partition induced by the signal θ∗, but also fulfill the target for (two-sided) consistency specified in (3).
  The two-step estimator. Our two-estimator starts with a constrained DCART, prohibiting splits resulting in rectangles that are too small. The first step estimation is defined as
  θ̂ = OS(Π̂)(y), with Π̂ ∈ arg min Π∈Pdyadic,d,n(η)
  { 2−1‖y −OS(Π)(y)‖2 + λ1|Π| } , (8)
  where λ1, η > 0 are tuning parameters and Pdyadic,d,n(η) is the set of rectangular partitions where every rectangle is of size at least η.
  The second step merges rectangles in the partition associated with the estimator from the first step to overcome over-partitioning. To be specific, let {Rl}l∈[k(θ̂)] be a rectangular partition of Ld,n
  associated by θ̂. For each (i, j) ∈ [k(θ̂)]× [k(θ̂)], i < j, let z(i,j) = 1 if
  dist(Ri, Rj) ≤ η (9)
  and 1
  2 ∑ l∈Ri (Yl − ȲRi)2 + ∑ l∈Rj (Yl − ȲRj )2 + λ2 > 1 2 ∑ l∈Ri∪Rj (Yl − ȲRi∪Rj )2, (10)
  where λ2 > 0 is a tuning parameter; otherwise, we let z(i,j) = 0. With this notation, letE = {(i, j) ∈ [k(θ̂)]× [k(θ̂)] : z(i,j) = 1} and let {Ĉl}l∈[L] be the collection of all the connected components of the undirected graph ([k(θ̂)], E). The final output can be written as
  Λ̂ = { ∪j∈Ĉ1Rj , . . . ,∪j∈ĈLRj } . (11)
  Notice that the main computational burden is to find the DCART estimator which has a cost of O(N). From the output of DCART, the computation of the quantities in (9) and (10) can be done in O(k(θ̂)2). Before describing the favorable properties of the estimator, we state our main assumption.
  Assumption 1 If A,B ∈ Λ∗ with A 6= B and θ̄∗A = θ̄∗B , then we have that
  dist(A,B) ≥ ckdyad(θ∗)κ−2σ2 log(N), (12)
  for some large enough constant c > 0. Furthermore, we assume that
  κ2∆ ≥ ckdyad(θ∗)σ2 log(N). (13)
  Assumption 1 simply states that if two elements of Λ∗ have the same signal values then they should be sufficiently apart from each other. Assumption 1 also specifies a signal-to-noise ratio type of condition. When d = 1, it is well known in the change point literature [e.g. 36] that the optimal signal-to-noise ratio for localization is of the form κ2∆ & σ2 log(N). The condition in (13) has an additional kdyad(θ∗) factor. It is an interesting open problem to determine whether this additional term can be avoided when d ≥ 2.
  Theorem 2 Assume that Assumption 1 holds. Suppose that the data satisfy (1) and Λ̂ is the two-step estimator, with tuning parameters λ1 = C1σ2 log(N), λ2 = C2kdyad(θ∗)σ2 log(N) and
  c1kdyad(θ ∗)σ2κ−2 log(N) ≤ η ≤ ∆/c2, (14)
  where c1, c2, C1, C2 > 0 are absolute constants. Then with probability at least 1−N−c, it holds that
  |Λ̂| = |Λ∗| and dHaus(Λ̂,Λ∗) ≤ Cσ2κ−2kdyad(θ∗) log(N), (15)
  where c, C > 0 are absolute constants.
  Theorem 2 shows that the two-step estimator overcomes the over-partitioning issue of DCART and is consistent for the partition recovery problem provided that κ−2σ2kdyad(θ∗) log(N)/∆ → 0. The resulting estimation error is of order at most κ−2σ2kdyad(θ∗) log(N).
  In view of Assumption 1 and Theorem 2, intuitively speaking, (14) ensures that if there are two separated regions where the true signal takes on the same value, then the they should be far apart; otherwise, our algorithm would not be able to tell if they should be merged together or keep separated. Eq. (15) requires that the signal strength is large enough. Technically speaking, if (14) is changed to another quantity, denoted by w, then the final result of Theorem 2 would be
  dHaus(Λ̂,Λ ∗) ≤ C { (w ∧∆) ∨ σ2κ−2kdyad(θ∗) log(N) } ,
  where the term w ∧∆ is due to the definition of ∆. This means that the final rate in Theorem 2 is determined jointly by Assumption 1 and an optimal rate.
  There are three tuning parameters required by the two-step estimator. Practical guidance on how to pick them will be provided in Section S6.2. The tuning parameter λ1 in (8) is the same as the one in (4) and their theoretical rates are determined by the maximal noise level of a Gaussian distribution over all possible rectangles using a union bound argument. The tuning parameter λ2 is used in
  the merging step (10), penalizing over-partitioning. Since the candidate rectangles in (10) are the estimators from the first step, these rectangles carry the estimation error from the first step. An intermediate result from the proof of Theorem 2 unveils a similar result to Theorem 1, that for each Rj involved in (10), there exists a subset Sj ⊂ Rj having constant signal values satisfying that∑
  j∈[k(θ̂)]
  |Rj\Sj | . κ−2σ2kdyad(θ∗) log(N).
  This suggests that the right choice for λ2 should be able to counter this extra kdyad(θ∗) factor. Finally, the tuning parameter η appears twice in the estimation procedure: as a lower bound on the size of the rectangles obtained in first step as (8) and as an upper bound on the distance between two rectangles in (9); see (14). The value of η should be at least as large as the one-sided upper bound on the recovery error, in order to ensure that over-partitioning cannot occur. As the same time, it should not be chosen too large, or otherwise the procedure may erroneously prune small true rectangles. By this logic, η should not exceed the minimal size of the true rectangles; this is indeed the upper bound in condition (14). Finally, we would like to point out that, similar conditions are even necessary in some change point localization (d = 1) procedures, see [e.g. 36].
  In practice, one may be tempted to abandon the tuning parameter η, and only prune the DCART output using (10). If one still wants the result to satisfy (3), then stronger conditions are needed and worse localization rates are obtained. We include this result in Section S4 in the supplementary material.
sections.heading: 2.3 Optimality in the regular boundary case
sections.text: We use a two-step procedure to improve the partition recovery performances of DCART and show the error is of order κ−2σ2kdyad(θ∗) log(N). A natural question in order is whether one can further expect to improve this rate. In Proposition 3, we show a minimax lower bound result.
  Proposition 3 Let {yi}i∈Ld,n satisfy (1) and
  θ∗i = κ, if i ∈ S; θ∗i = 0, if i ∈ Ld,n \ S, (16)
  where S ⊂ Ld,n is a rectangle and |S| = ∆ > 0. Let PNκ,∆,σ denote the corresponding joint distribution. Consider the class of distributions
  PN = { PNκ,∆,σ : ∆ < N/2, κ 2∆/σ2 ≥ log(N)/6 } .
  Then forN ≥ 26, it holds that inf Ŝ supP∈PN EP { |Ŝ4S| } ≥ σ2κ−2 log(N)/36, where the infimum
  is over all estimators Ŝ of S.
  Theorem 3 shows that when the induced partition of θ∗ consists of one rectangle and its complement, i.e. when kdyad(θ∗) = O(1), the minimax lower bound on the estimation error is of order κ−2σ2 log(N). Recalling the estimation errors we derived for DCART and the two-step estimator in Theorems 1 and 2, when kdyad(θ∗) = O(1), the results thereof are minimax optimal.
  The assumption kdyad(θ∗) = O(1) is fairly restrictive, though, using our notation, the case of |Λ∗| = 2 is in fact used in the cluster detection literature [4, 3]. In fact, in order to achieve the optimal estimation rate indicated in Theorem 3, we only need a boundary regularity condition, in the sense that for every rectangle in the partition induced by θ∗, there are only O(1)-many other rectangles nearby. This condition is formally stated in Assumption 2.
  Assumption 2 There exists constant C, c > 0 such that for any A ∈ Λ∗ it holds that∣∣{B ∈ Λ∗\{A} : dist(A,B) ≤ cσ2κ−2kdyad(θ∗) log(N) }∣∣ ≤ C. Assumption 2 asserts that within cκ−2σ2kdyad(θ∗) log(N) distance, each element of Λ∗ only has O(1)-many elements nearby. This condition shares the same spirit of requiring the cluster boundary to be a bi-Lipschitz function in the cluster detection literature [e.g. 2].
  Corollary 4 Assume that Assumptions 1 and 2 hold. Suppose that the data satisfy (1) and Λ̂ is the two-step estimator defined in (11), with tuning parameters λ1 = C1σ2 log(N), λ2 = C2kdyad(θ
  ∗)σ2 log(N) and c1κ−2kdyad(θ∗)σ2 log(N) ≤ η ≤ ∆/c2, where c1, c2, C1, C2 > 0 are absolute constants. Then with probability at least 1−N−c, it holds that |Λ̂| = |Λ∗| and
  max A∈Λ∗ min Â∈Λ̂ |Â4A∗| ≤ Cσ
  2 log(N)
  κ2 min
  { k(θ∗), max
  B∈Λ∗, dist(A,B)≤cσ2κ−2kdyad(θ∗) log(N) |B|/η
  } .
  (17) where c, C > 0 are absolute constants. In particular, if η ∆ and |A| ∆, for all A ∈ Λ∗, then
  dHaus(Λ̂,Λ ∗) ≤ C σ2κ−2 log(N). (18)
  Corollary 4 shows that even if kdyad(θ∗) is diverging as the sample size grows unbounded, one can still achieve the minimax optimal estimation error rate κ−2σ2 log(N), with properly chosen tuning parameters and additional regularity conditions on the partition. An interesting by-product in deriving this rate is (17), which characterizes the effect of the number of nearby rectangles in the estimation error for individual elements in Λ∗.
sections.heading: 3 Experiments
sections.text: In this section, we demonstrate in simulation the numerical performances of the two-step estimator for the task of partition recovery. The code is by courtesy of the authors of [12] and all of our experiments are done in a 2.3 GHz 8-Core Intel Core i9 machine. Our code can be found in https://github. com/hernanmp/Partition_recovery. We focus on the naive two-step estimator detailed in Section S4 and denoted here as Λ̂. The implementation details regarding choice of tuning parameters are discussed in S6.2.
  We adopt dist1 = dHaus(Λ̂,Λ∗) and dist2 = ||Λ̂| − |Λ∗|| as the measurements. For each scenario depicted in Figure 3, we report the mean and standard errors of dist1 and dist2 over 50 Monte Carlo simulations.
  As a competitor benchmark, we consider a similar pruning algorithm based on the total variation estimator [27, 33], namely TV-based, instead of ours based on DCART. The implementation details are discussed in Section S6.3.
  For each scenario considered, we vary the noise level as σ ∈ {0.5, 1, 1.5} and set (d, n) = (2, 27). In each instance, the data are generated as y ∼ N (θ∗, σ2ILd,n). Detailed descriptions are in Section S6.1, and visualizations of the signal patterns are in the second column in Figure 3, while the third and fourth columns depict θ̃, the DCART estimator, and θ̂, our two-step estimator, respectively. We can see that our two-step estimator correctly identifies the partition and improves upon DCART for the purpose of recovery partition. It is worth noting that even when the partition is not rectangular, as shown in the second row in Figure 3, our two-step estimator is still able to accurately recover a good rectangular partition.
  From Table 1 we see that in terms of the metric dist2 our two-step estimator outperforms TV-based estimator in all cases. Furthermore, the same is also true for most cases in terms of the metric dist1.
sections.heading: 4 Conclusions
sections.text: In this paper we study the partition recovery problem over d-dimensional lattices. We show how a simple modification of DCART enjoys one-sided consistency. To correctly identify the size of the true partition and obtain better consistency guarantees, we further propose a more sophisticated two-step estimator which is shown to be minimax optimal in some regular cases.
  Throughout this paper, we discuss partition recovery on lattice grids. In fact, to deal with non axis-aligned data, one can construct a lattice in the domain of the features, average the observations within each cell of the lattice and ignore the cells without observations. More details can be found in the supplementary materials.
  Finally, an important open question remains regarding the necessity to include the factor kdyad(θ∗) in the signal-to-noise condition (13) and in the estimation rates. We leave for future work to investigate what the optimal estimation rates are in general when kdyad(θ∗) is allowed to diverge.
  Supplementary material
  The supplementary material contains all the technical details, additional simulation results and code of this paper.
  Acknowledgments and Disclosure of Funding
  Funding in direct support of this work: NSF DMS 2015489 and EPSRC EP/V013432/1.
title: Lattice partition recovery with dyadic CART