summaries: This paper proposes a way of learning a capsule-based representation of 3D point clouds. The representation decomposes an input point cloud into parts, each corresponding to a capsule. This decomposition is trained to be consistent across SE(3) transformations of the point cloud using loss terms that encourage feature invariance and pose equivariance across two random presentations of the input point cloud. Moreover, the model learns a canonical pose for each part that is tied to the identity of the object (not its presented pose).  The experiments show that 1. The reconstruction quality is high, showing that the model has capacity to fit point clouds from ShapeNet. 1. The canonical pose is indeed consistent (measured by low std deviation in relative angle of rotation). 1. The model can be used for pairwise registration (but DeepGMR does even better using RRI features, but doesn't do any of the other things this model can do.) 1. The learned features lend themselves well to unsupervised clustering/classification. 1. All the proposed loss terms are useful and contribute to better reconstruction and canonicalization.

summaries: The paper proposes a capsule-based network architecture that allows self-supervised representation learning from 3D point clouds. While existing methods rely heavily on pre-canonicalized (aligned) 3D models for training, this work proposes to automatically discover the "optimal" canonical transformation for alignment while learning to decompose and reconstruct the input point clouds.  Specifically, it exploits the transformation invariance and equivariance properties of object parts (capsules) to learn a latent representation. The experimental results on the (unaligned) ShapeNet dataset show that the method can effectively learn to decompose the 3D point cloud into consistent parts across objects and canonicalize unaligned data. It outperforms the state-of-the-art methods in terms of reconstruction accuracy, canonicalization accuracy and stability, as well as unsupervised classification accuracy from latent representation. It is very encouraging to see that the need of 3D pose/viewpoint annotations can be removed and such competitive performance is reached.

summaries: This paper proposes to learn a Canonical Capsule decomposition of shape point clouds in a self-supervised manner.  Given an unaligned 3D point cloud dataset, the proposed method constructs pairwise data with random SE(3) transformations, and then extracts the per-capsule pose and a transformation-invariant descriptor. The experiment results show the importance of the object-centric representation with capsule decomposition in many applications, including reconstruction, canonicalization, and classification.

summaries: The paper proposes a self-supervised capsule architecture for 3D point clouds. In order to circumvent the customary need to train on pre-aligned datasets, the authors introduce the Canonical Capsules to compute the K-part decomposition of a point cloud that allows unsupervised "object-centric" learning of 3D representations without requiring a semantically aligned dataset. The main idea of the canonical capsules is to enforce invariance/equivariance by relating the decomposition result of two random rigid transformations.  The authors compare the proposed method against state-of-the-art methods on 3D point cloud reconstruction, canonicalization, and unsupervised classification.

summaries: The paper proposed a self-supervised capsule based architecture for 3D point cloud. The main goal is to learn a  canonicalization operation that allows unbiased object-centric reasoning.  The idea is that an object is decomposed to K keypoints using canonical capsules, the whole system is trained in an unsupervised fashion based on pairs of randomly rotated/translated  of the same object. The capsule pose equivariance is enforced by requiring the equivariance of two keypoint sets based on the known relative transformation of the pair images and the capsule descriptor invariance is enforced by matching the pair keypoints. The proposed framework have been evaluated in 3D autoencoding and reconstruction, unsupervised classification and canonicalization on ShapeNet dataset.

abstractText: We propose a self-supervised capsule architecture for 3D point clouds. We compute capsule decompositions of objects through permutation-equivariant attention, and self-supervise the process by training with pairs of randomly rotated objects. Our key idea is to aggregate the attention masks into semantic keypoints, and use these to supervise a decomposition that satisfies the capsule invariance/equivariance properties. This not only enables the training of a semantically consistent decomposition, but also allows us to learn a canonicalization operation that enables object-centric reasoning. To train our neural network we require neither classification labels nor manually-aligned training datasets. Yet, by learning an object-centric representation in a self-supervised manner, our method outperforms the state-of-the-art on 3D point cloud reconstruction, canonicalization, and unsupervised classification.

title: Pattern Recognition and Machine Learning

venue: springer,

authors:
Rohan Chabra
Jan Eric Lenssen
Eddy Ilg
Tanner Schmidt
Julian Straub
Steven Lovegrove
Richard Newcombe

title: Deep Local Shapes: Learning Local SDF Priors for Detailed 3D Reconstruction

venue: In European Conference on Computer Vision,

authors:
Angel X Chang
Thomas Funkhouser
Leonidas Guibas
Pat Hanrahan
Qixing Huang
Zimo Li
Silvio Savarese
Manolis Savva
Shuran Song
Hao Su

title: Shapenet: An Information-Rich 3D Model Repository

venue: arXiv Preprint,

authors:
Chao Chen
Guanbin Li
Ruijia Xu
Tianshui Chen
Meng Wang
Liang Lin

title: Clusternet: Deep Hierarchical Cluster Network with Rigorously Rotation-Invariant Representation for Point Cloud Analysis

venue: In Conference on Computer Vision and Pattern Recognition,

authors:
Ting Chen
Simon Kornblith
Mohammad Norouzi
Geoffrey Hinton

title: A Simple Framework for Contrastive Learning of Visual Representations

venue: International Conference on Machine Learning,

authors:
Zhiqin Chen
Andrea Tagliasacchi
Hao Zhang

title: Bsp-net: Generating compact meshes via binary space partitioning

authors:
Zhiqin Chen
Hao Zhang

title: Learning Implicit Fields for Generative Shape Modeling

authors:
Boyang Deng
JP Lewis
Timothy Jeruzalski
Gerard Pons-Moll
Geoffrey Hinton
Mohammad Norouzi
Andrea Tagliasacchi

title: NASA: Neural Articulated Shape Approximation

authors:
Congyue Deng
Or Litany
Yueqi Duan
Adrien Poulenard
Andrea Tagliasacchi
Leonidas Guibas

title: Vector neurons: a general framework for so(3)-equivariant networks, 2021

authors:
Zhantao Deng
Jan Bednařík
Mathieu Salzmann
Pascal Fua

title: Better Patch Stitching for Parametric Surface Reconstruction

authors: Deng Boyang Genova Kyle Yazdani Soroosh Bouaziz Sofien Hinton Geoffrey Tagliasacchi Andrea

title: CvxNet: Learnable Convex Decomposition

authors:
Theo Deprelle
Thibault Groueix
Matthew Fisher
Vladimir Kim
Bryan Russell
Mathieu Aubry

title: Learning Elementary Structures for 3D Shape Generation and Matching

venue: In Advances in Neural Information Processing Systems,

authors:
Haoqiang Fan
Hao Su
Leonidas J Guibas

title: A Point Set Generation Network for 3D Object Reconstruction from a Single Image

authors:
Pedro Felzenszwalb
David McAllester
Deva Ramanan

title: A Discriminatively Trained, Multiscale

venue: Deformable Part Model. In Conference on Computer Vision and Pattern Recognition,

authors:
Clara Fernandez-Labrador
Ajad Chhatkuli
Danda Pani Paudel
Jose J Guerrero
Cédric Demonceaux
Luc Van Gool

title: Unsupervised Learning of Category-Specific Symmetric 3D Keypoints from Point Sets

authors:
Kyle Genova
Forrester Cole
Avneesh Sud
Aaron Sarna
Thomas Funkhouser

title: Deep Structured Implicit Functions

authors:
Kyle Genova
Forrester Cole
Daniel Vlasic
Aaron Sarna
William T Freeman
Thomas Funkhouser

title: Learning Shape Templates with Structured Implicit Functions

venue: In International Conference on Computer Vision,

authors:
Zan Gojcic
Caifa Zhou
Jan D Wegner
Leonidas J Guibas
Tolga Birdal

title: Learning multiview 3D point cloud registration

authors:
Thibault Groueix
Matthew Fisher
Vladimir G Kim
Bryan C Russell
Mathieu Aubry

title: A Papier-Mâché Approach to Learning 3D Surface Generation

authors:
Jiayuan Gu
Wei-Chiu Ma
Sivabalan Manivasagam
Wenyuan Zeng
Zihao Wang
Yuwen Xiong
Hao Su
Raquel Urtasun

title: Weakly-Supervised 3D Shape Completion in the Wild

authors:
Kaiming He
Xiangyu Zhang
Shaoqing Ren
Jian Sun

title: Deep Residual Learning for Image Recognition

authors:
Geoffrey E Hinton
Alex Krizhevsky
Sida D Wang

title: Transforming Auto-Encoders

venue: In International Conference on Artificial Neural Networks,

authors:
Geoffrey E Hinton
Kevin J Lang

title: Shape Recognition and Illusory Conjunctions

venue: In International Joint Conference on Artificial Intelligence,

title: A Parallel Computation that Assigns Canonical Object-based Frames of Reference

authors:
Sergey Ioffe
Christian Szegedy

title: Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift

venue: In International Conference on Machine Learning,

authors:
Arthur Jacot
Franck Gabriel
Clément Hongler

title: Neural Tangent Kernel: Convergence and Generalization in Neural Networks

authors:
Abhishek Kar
Christian Häne
Jitendra Malik

title: Learning a Multi-View Stereo Machine

authors:
D.P. Kingma
J. Ba

title: Adam: A Method for Stochastic Optimisation

venue: In International Conference on Learning Representations,

authors:
Adam Kosiorek
Sara Sabour
Yee Whye Teh
Geoffrey E Hinton

title: Stacked Capsule Autoencoders

authors:
Or Litany
Alex Bronstein
Michael Bronstein
Ameesh Makadia

title: Deformable Shape Completion with Graph Convolutional Autoencoders

title: Distinctive Image Features from Scale-Invariant Keypoints

venue: International Journal of Computer Vision,

authors:
Lars Mescheder
Michael Oechsle
Michael Niemeyer
Sebastian Nowozin
Andreas Geiger

title: Occupancy Networks: Learning 3D Reconstruction in Function Space

authors:
David Novotny
Nikhila Ravi
Benjamin Graham
Natalia Neverova
Andrea Vedaldi

title: C3dpo: Canonical 3d pose networks for non-rigid structure from motion

authors:
Jeong Joon Park
Peter Florence
Julian Straub
Richard Newcombe
Steven Lovegrove

title: DeepSDF: Learning Continuous Signed Distance Functions for Shape Representation

authors:
Charles R Qi
Hao Su
Kaichun Mo
Leonidas J Guibas

title: Pointnet: Deep Learning on Point Sets for 3D Classification and Segmentation

authors:
Charles R Qi
Li Yi
Hao Su
Leonidas J Guibas

title: PointNet++: Deep Hierarchical Feature Learning on Point Sets in a Metric Space

venue: Advances in Neural Information Processing Systems,

authors:
Joseph Redmon
Santosh Divvala
Ross Girshick
Ali Farhadi

title: You Only Look Once: Unified, Real-Time Object Detection

authors:
Davis Rempe
Tolga Birdal
Yongheng Zhao
Zan Gojcic
Srinath Sridhar
Leonidas J. Guibas

title: CaSPR: Learning Canonical Spatiotemporal Point Cloud Representations

authors:
Sara Sabour
Andrea Tagliasacchi
Soroosh Yazdani
Geoffrey E. Hinton
David J. Fleet

title: Unsupervised part representation by flow capsules, 2021

title: Uniform random rotations

venue: In Graphics Gems III (IBM Version),

authors:
Karen Simonyan
Andrew Zisserman

title: Very Deep Convolutional Networks for Large-Scale Image Recognition

venue: International Conference on Learning Representations,

authors:
Shuran Song
Fisher Yu
Andy Zeng
Angel X Chang
Manolis Savva
Thomas Funkhouser

title: Semantic scene completion from a single depth image

venue: Proceedings of 30th IEEE Conference on Computer Vision and Pattern Recognition,

authors:
Olga Sorkine-Hornung
Michael Rabinovich

title: Least-Squares Rigid Motion

venue: Using SVD. Computing,

authors:
Riccardo Spezialetti
Federico Stella
Marlon Marcon
Luciano Silva
Samuele Salti
Luigi Di Stefano

title: Learning to Orient Surfaces by Self-supervised Spherical CNNs

authors:
Nitish Srivastava
Hanlin Goh
Ruslan Salakhutdinov

title: Geometric Capsule Autoencoders for 3D Point Clouds

authors:
Weiwei Sun
Wei Jiang
Eduard Trulls
Andrea Tagliasacchi
Kwang Moo Yi

title: ACNe: Attentive Context Normalization for Robust Permutation-Equivariant Learning

authors:
Supasorn Suwajanakorn
Noah Snavely
Jonathan J Tompson
Mohammad Norouzi

title: Discovery of Latent 3D Keypoints via End-to-End Geometric Reasoning

authors:
Nathaniel Thomas
Tess Smidt
Steven Kearnes
Lusann Yang
Li Li
Kai Kohlhoff
Patrick Riley

title: Tensor Field Networks: Rotation-and Translation-Equivariant Neural Networks for 3D Point Clouds

authors:
Shubham Tulsiani
Hao Su
Leonidas J Guibas
Alexei A Efros
Jitendra Malik

title: Learning shape abstractions by assembling volumetric primitives

authors:
Oliver van Kaick
Kai Xu
Hao Zhang
Yanzhen Wang
Shuyang Sun
Ariel Shamir
Daniel Cohen-Or

title: Co-hierarchical analysis of shape structures

venue: ACM SIGGRAPH,

authors:
He Wang
Srinath Sridhar
Jingwei Huang
Julien Valentin
Shuran Song
Leonidas J Guibas

title: Normalized Object Coordinate Space for Category-level 6d Object Pose and Size Estimation

authors:
Nanyang Wang
Yinda Zhang
Zhuwen Li
Yanwei Fu
Wei Liu
Yu-Gang Jiang

title: Pixel2mesh: Generating 3d mesh models from single rgb images

authors:
Peng-Shuai Wang
Yang Liu
Yu-Xiao Guo
Chun-Yu Sun
Xin Tong

title: O-CNN: Octree-Based Convolutional Neural Networks for 3d Shape Analysis

venue: ACM Transactions on Graphics,

authors:
Peng-Shuai Wang
Chun-Yu Sun
Yang Liu
Xin Tong

title: Adaptive o-cnn: A patch-based deep representation of 3d shapes

venue: ACM Transactions on Graphics (TOG),

authors:
Yue Wang
Justin M Solomon

title: Deep Closest Point: Learning Representations for Point Cloud Registration

authors:
Yue Wang
Yongbin Sun
Ziwei Liu
Sanjay E Sarma
Michael M Bronstein
Justin M Solomon

title: Dynamic graph cnn for learning on point clouds

authors:
Yuxin Wu
Kaiming He

title: Group Normalization

authors:
Saining Xie
Jiatao Gu
Demi Guo
Charles R. Qi
Leonidas J. Guibas
Or Litany

title: PointContrast: Unsupervised Pre-training for 3D Point Cloud Understanding

authors:
Xinchen Yan
Jimei Yang
Ersin Yumer
Yijie Guo
Honglak Lee

title: Perspective Transformer Nets: Learning Single-View 3D Object Reconstruction without 3D Supervision

authors:
Yaoqing Yang
Chen Feng
Yiru Shen
Dong Tian

title: FoldingNet: Point Cloud Auto-Encoder via Deep Grid Deformation

authors:
Zi Jian Yew
Gim Hee Lee

title: Rpm-net: Robust point matching using learned features

authors:
Wentao Yuan
Ben Eckart
Kihwan Kim
Varun Jampani
Dieter Fox
Jan Kautz

title: DeepGMR: Learning Latent Gaussian Mixture Models for Registration

authors:
Yongheng Zhao
Tolga Birdal
Haowen Deng
Federico Tombari

title: 3D Point Capsule Networks

authors:
Yongheng Zhao
Tolga Birdal
Jan Eric Lenssen
Emanuele Menegatti
Leonidas Guibas
Federico Tombari

title: Quaternion Equivariant Capsule Networks for 3D Point Clouds

heading: 1 Introduction

text:
Understanding objects is one of the core problems of computer vision [32, 14, 38]. While this task has traditionally relied on large annotated datasets [42, 22], unsupervised approaches that utilize self-supervision [5] have emerged to remove the need for labels. Recently, researchers have attempted to extend these methods to work on 3D point clouds [59], but the field of unsupervised 3D learning remains relatively uncharted. Conversely, researchers have been extensively investigating 3D deep representations for shape auto-encoding1 [61, 19, 33, 16], making one wonder whether these discoveries can now benefit from unsupervised learning for tasks other than auto-encoding.
Importantly, these recent methods for 3D deep representation learning are not entirely unsupervised. Whether using point clouds [61], meshes [19], or implicits [33], they owe much of their success to the bias within the dataset that was used for training. Specifically, all 3D models in the popular ShapeNet [3] dataset are “object-centric” – they are pre-canonicalized to a unit bounding box, and, even more importantly, with an orientation that synchronizes object semantics to Euclidean frame axes (e.g. airplane cockpit is always along +y, car wheels always touch z = 0). Differentiable 3D decoders are heavily affected by the consistent alignment of their output with an Euclidean frame [8, 16] as local-to-global transformations cannot be easily learnt by fully connected layers. As we will show in Section 4.2, these methods fail in the absence of pre-alignment, even when data
1Auto-encoding is also at times referred to as “reconstruction” or “shape-space” learning.
35th Conference on Neural Information Processing Systems (NeurIPS 2021), virtual.
augmentation is used. A concurrent work [45] also recognizes this problem and proposes a separate learnt canonicalizer, which is shown helpful in downstream classification tasks.
In this work, we leverage the modeling paradigm of capsule networks [23]. In capsule networks, a scene is perceived via its decomposition into part hierarchies, and each part is represented with a (pose, descriptor) pair: 1 The capsule pose specifies the frame of reference of a part, and hence should be transformation equivariant; 2 The capsule descriptor specifies the appearance of a part, and hence should be transformation invariant. Thus, one does not have to worry how the data is oriented or translated, as these changes can be encoded within the capsule representation.
We introduce Canonical Capsules, a novel capsule architecture to compute a K-part decomposition of a point cloud. We train our network by feeding pairs of a randomly rotated/translated copies of the same shape (i.e. siamese training) hence removing the requirement of pre-aligned training datasets. We then decompose the point cloud by assigning each point into one of the K parts via attention, which we aggregate into K keypoints. Equivariance is then enforced by requiring the two keypoint sets to only differ by the known (relative) transformation (i.e. a form of self-supervision). For invariance, we simply ask that the descriptors of each keypoint of the two instances match.
With Canonical Capsules, we exploit our decomposition to recover a canonical frame that allows unsupervised “object-centric” learning of 3D deep representations without requiring a semantically aligned dataset. We achieve this task by regressing canonical capsule poses from capsule descriptors via a deep network, and computing a canonicalizing transformation by solving a shape-matching problem [44]. This not only allows more effective shape auto-encoding, but our experiments confirm this results in a latent representation that is more effective in unsupervised classification tasks. Note that, like our decomposition, our canonicalizing transformations are also learnt in a self-supervised fashion, by only training on randomly transformed point clouds.
Contributions. In summary, in this paper we:
• propose an architecture for 3D self-supervised learning based on capsule decomposition; • enable object-centric unsupervised learning by introducing a learned canonical frame of reference; • achieve state-of-the-art performance 3D point cloud auto-encoding/reconstruction, canonicalization,
and unsupervised classification.

heading: 2 Related works

text:
Convolutional Neural Networks lack equivariance to rigid transformations, despite their pivotal role in describing the structure of the 3D scene behind a 2D image. One promising approach to overcome this shortcoming is to add equivariance under a group action in each layer [49, 9]. In our work, we remove the need for a global SE(3)-equivariant network by canonicalizing the input.
Capsule Networks. Capsule Networks [23] have been proposed to overcome this issue towards a relational and hierarchical understanding of natural images. Of particular relevance to our work, are methods that apply capsule networks to 3D input data [64, 65, 46], but note these methods are not unsupervised, as they either rely on classification supervision [65], or on datasets that present a significant inductive bias in the form of pre-alignment [64]. In this paper, we take inspiration from the recent Stacked Capsule Auto-Encoders [30], which shows how capsule-style reasoning can be effective as long as the primary capsules can be trained in a self-supervised fashion (i.e. via reconstruction losses). The natural question, which we answer in this paper, is “How can we engineer networks that generate 3D primary capsules in an unsupervised fashion?”
Deep 3D representations. Reconstructing 3D objects requires effective inductive biases about 3D vision and 3D geometry. When the input is images, the core challenge is how to encode 3D projective geometry concepts into the model. This can be achieved by explicitly modeling multi-view geometry [28], by attempting to learn it [13], or by hybrid solutions [60]. But even when input is 3D, there are still significant challenges. It is still not clear which is the 3D representation that is most amenable to deep learning. Researchers proposed the use of meshes [53, 31], voxels [54, 55], surface patches [19, 12, 10], and implicit functions [35, 33, 7]. Unfortunately, the importance of geometric structures (i.e. part-to-whole relationships) is often overlooked. Recent works have tried to close this gap by using part decomposition consisting of oriented boxes [50], ellipsoids [17, 16], convex polytopes [11], and grids [2]. However, as previously discussed, most of these still heavily rely on a
pre-aligned training dataset; our paper attempts to bridge this gap, allowing learning of structured 3D representations without requiring pre-aligned data.
Canonicalization. One way to circumvent the requirement of pre-aligned datasets is to rely on methods capable of registering a point cloud into a canonical frame. The recently proposed CaSPR [39] fulfills this premise, but requires ground-truth canonical point clouds in the form of normalized object coordinate spaces [52] for supervision. Similarly, [20] regresses each view’s pose relative to the canonical pose, but still requires weak annotations in the form of multiple partial views. C3DPO [34] learns a canonical 3D frame based on 2D input keypoints. In contrast to these methods, our solution is completely self-supervised. The concurrent Compass [45] also learns to canonicalize in a self-supervised fashion, but as the process is not end-to-end, this results in a worse performance than ours, as it will be shown in Section 4.3.
Registration. Besides canonicalization, many pairwise registration techniques based on deep learning have been proposed (e.g. [56, 63]), even using semantic keypoints and symmetry to perform the task [15]. These methods typically register a pair of instances from the same class, but lack the ability to jointly and consistently register all instances to a shared canonical frame.

heading: 3 Method

text:
Our network trains on unaligned point clouds as illustrated in Figure 1: we train a network that decomposes point clouds into parts, and enforce invariance/equivariance through a Siamese training setup [48]. We then canonicalize the point cloud to a learnt frame of reference, and perform autoencoding in this coordinate space. The losses employed to train E , K, and D, will be covered in Section 3.1, while the details of their architecture are in Section 3.2.
Decomposition. In more detail, given a point cloud P 2 RP⇥D of P points in D dimensions (in our case D=3), we perturb it with two random transformations Ta,Tb 2 SE(D) to produce point clouds Pa,Pb. We then use a shared permutation-equivariant capsule encoder E to compute a K-fold attention map A 2 RP⇥K for K capsules, as well as per-point feature map F 2 RP⇥C with C channels: A,F = E(P) , (1) where we drop the superscript indexing the Siamese branch for simplicity. From these attention masks, we then compute, for the k–th capsule its pose ✓k 2 R3 parameterized by its location in 3D space, and the corresponding capsule descriptor k 2 RC :
✓k = P p Ap,kPpP p Ap,k , k = P p Ap,kFpP p Ap,k . (2)
Hence, as long as E is invariant w.r.t. rigid transformations of P, the pose ✓k will be transformation equivariant, and the descriptor k will be transformation invariant. Note that this simplifies the design (and training) of the encoder E , which only needs to be invariant, rather than equivariant [49, 46]. Canonicalization. Simply enforcing invariance and equivariance with the above framework is not enough to learn 3D representations that are object-centric, as we lack an (unsupervised) mechanism to bring information into a shared “object-centric” reference frame. Furthermore, the “right” canonical
frame is nothing but a convention, thus we need a mechanism that allows the network to make a choice – a choice, however, that must then be consistent across all objects. For example, a learnt canonical frame where the cockpit of airplanes is consistently positioned along +z is just as good as a canonical frame where it is positioned along the +y axis. To address this, we propose to link the capsule descriptors to the capsule poses in canonical space, that is, we ask that objects with similar appearance to be located in similar Euclidean neighborhoods in canonical space. We achieve this by regressing canonical capsules poses (i.e. canonical keypoints) ✓̄ 2 RK⇥3 using the descriptors 2 RK⇥C via a fully connected deep network K:
✓̄ = K ( ) (3) Because fully connected layers are biased towards learning low-frequency representations [27], this regressor also acts as a regularizer that enforces semantic locality.
Auto-encoding. Finally, in the learnt canonical frame of reference, to train the capsule descriptors via auto-encoding, we reconstruct the point clouds with per-capsule decoders Dk:
P̃ = [k Dk(R̄✓k + t̄, k) , (4)
where [ denotes the union operator. The canonicalizing transformation T̄ = (R̄, t̄) can be readily computed by solving a shape-matching problem [44], thanks to the property that our capsule poses and regressed keypoints are in one-to-one correspondence:
R̄, t̄ = argmin R,t
1
K
X
k
k(R✓k + t) ✓̄kk22 . (5)
While the reconstruction in (4) is in canonical frame, note it is trivial to transform the point cloud back to the original coordinate system after reconstruction, as the transformation T̄ 1 is available.

heading: 3.1 Losses

text:
As common in unsupervised methods, our framework relies on a number of losses that control the different characteristics we seek to obtain in our representation. Note none of these losses require ground truth labels. We organize the losses according to the portion of the network they supervise: decomposition, canonicalization, and reconstruction.
Decomposition. While a transformation invariant encoder architecture should be sufficient to achieve equivariance/invariance, this does not prevent the encoder from producing trivial solutions/decompositions once trained. As capsule poses should be transformation equivariant, the poses of the two rotation augmentations ✓ak and ✓ b k should only differ by the (known) relative transformation:
Lequivariance = 1
K
X
k
k✓ak (Ta)(Tb) 1✓ b kk22 . (6)
Conversely, capsule descriptors should be transformation invariant, and as the two input points clouds are of the same object, the corresponding capsule descriptors should be identical:
Linvariance = 1
K
X
k
k ak b kk22 . (7)
We further regularize the capsule decomposition to ensure each of the K heads roughly represent the same “amount” of the input point cloud, hence preventing degenerate (zero attention) capsules. This is achieved by penalizing the attention variance:
Lequilibrium = 1
K
X
k
kak 1K⌃kakk 2 2 , (8)
where ak = ⌃p(Ap,k) denotes the total attention exerted by the k-th head on the point cloud.
Finally, to facilitate the training process, we ask for capsules to learn a localized representation of geometry. We express the spatial extent of a capsule by computing first-order moments of the represented points with respect to the capsule pose ✓k:
Llocalization = 1
K
X
k
1 ak
X
p
Ap,kk✓k Ppk22 . (9)
Canonicalization. To train our canonicalizer K, we relate the predicted capsule poses to regressed canonical capsule poses via the optimal rigid transformation from (5):
Lcanonical = 1
K
X
k
k(R̄✓k + t̄) ✓̄kk22 . (10)
Recall that R̄ and T̄ are obtained through a differentiable process. Thus, this loss is forcing the aggregated pose ✓k to agree with the one that goes through the regression path, ✓̄k. Now, since ✓̄k is regressed solely from the set of capsule descriptors, similar shapes will result in similar canonical keypoints, and the coordinate system of ✓̄k is one that employs Euclidean space to encode semantics.
Reconstruction. To learn canonical capsule descriptors in an unsupervised fashion, we rely on an auto-encoding task. We train the decoders {Dk} by minimizing the Chamfer Distance (CD) between the (canonicalized) input point cloud and the reconstructed one, as in [61, 19]:
Lrecon = CD ⇣ R̄P+ t̄, P̃ ⌘ . (11)

heading: 3.2 Network Architectures

text:
We briefly summarize our implementation details, including the network architecture; for further details, please refer to the supplementary material.
Encoder – E . Our architecture is based on the one suggested in [47]: a pointnet-like architecture with residual connections and attentive context normalization. We utilize Batch Normalization [26] instead of the Group Normalization [58], which trained faster in our experiments. We further extend their method to have multiple attention maps, where each attention map corresponds to a capsule.
Decoder – D. The decoder from (4) operates on a per-capsule basis. Our decoder architecture is similar to AtlasNetV2 [12] (with trainable grids). The difference is that we translate the per-capsule decoded point cloud by the corresponding capsule pose.
Regressor – K. We concatenate the descriptors and apply a series of fully connected layers with ReLU activation to regress the P capsule locations. At the output layer we use a linear activation and subtract the mean of the outputs to make our regressed locations zero-centered in the canonical frame.
Canonicalizing the descriptors. As our descriptors are only approximately rotation invariant (via Linvariance), we found it useful to re-extract the capsule descriptors k after canonicalization. Specifically, we compute F̄ with the same encoder setup, but with P̄=R̄P+T̄ instead of P and use it to compute ̄k; we validate this empirically in the supplementary material.

heading: 4 Results

text: We first discuss the experimental setup, and then validate our method on a variety of tasks: autoencoding, canonicalization, and unsupervised classification. While the task differs, our learning process remains the same: we learn capsules by reconstructing objects in a learnt canonical frame. We also provide an ablation study, which is expanded in detail in the supplementary material, where we provide additional qualitative results.

heading: 4.1 Experimental setup

text:
To evaluate our method, we rely on the ShapeNet (Core) dataset [3]2. We follow the category choices from AtlasNetV2 [12], using the airplane and chair classes for single-category experiments, while for multi-category experiments we use all 13 classes: airplane, bench, cabinet, car, chair, monitor, lamp, speaker, firearm, couch, table, cellphone, and watercraft. To make our results most compatible with those reported in the literature, we also use the same splits as in AtlasNetV2 [12]: 31747 shapes in the train, and 7943 shapes in the test set. Unless otherwise noted, we randomly sample 1024 points from the object surface for each shape to create our 3D point clouds.
De-canonicalizing the dataset. As discussed in the introduction, ShapeNet (Core) contains substantial inductive bias in the form of consistent semantic alignment. To remove this bias, we create random
2Available to researchers for non-commercial research and educational use.
SE(3) transformations, and apply them to each point cloud. We first generate uniformly sampled random rotations, and add uniformly sampled random translations within the range [ 0.2, 0.2], where the bounding volume of the shape ranges in [ 1,+1]. Note the relatively limited translation range is chosen to give state-of-the-art methods a chance to compete with our solution. We then use the relative transformation between the point clouds extracted from this ground-truth transformation to evaluate our methods. We refer to this unaligned version of the ShapeNet Core dataset as the unaligned setup, and using the vanilla ShapeNet Core dataset as the aligned setup. For the aligned setup, as there is no need for equivariance adaptation, we simply train our method without the random transformations, and so Lequivariance and Linvariance are not used. This setup is to simply demonstrate how Canonical Capsules would perform in the presence of a dataset bias.
We emphasize here that a proper generation of random rotation is important. While some existing works have generated them by uniformly sampling the degrees of freedom of an Euler-angle representation, this is known to be an incorrect way to sample random rotations [41], leading to biases in the generated dataset; see the supplementary material.
Implementation details. For all our experiments we use the Adam optimizer [29] with an initial learning rate of 0.001 and decay rate of 0.1. We train for 325 epochs for the aligned setup to match the AtlasNetV2 [12] original setup. For the unaligned setting, as the problem is harder, we train for a longer number of 450 epochs. We use a batch size of 16. The training rate is ⇠2.5 iters/sec. We train each model on a single NVidia V100 GPU. Unless stated otherwise, we use k=10 capsules and capsule descriptors of dimension C=128. We train three models with our method: two that are single-category (i.e., for airplane and chairs), and one that is multi-category (i.e., all 13 classes). To set the weights for each loss term, we rely on the reconstruction performance (CD) in the training set. We set weights to be one for all terms except for Lequivariance (5) and Lequilibrium (10 3). In the aligned case, because Lequivariance and Linvariance are not needed (always zero), we reduce the weights for the other decomposition losses by 103; Llocalization to 10 3 and Lequilibrium to 10 6.

heading: 4.2 Auto-encoding – Figure 2 and Table 1

text:
We evaluate the performance of our method for the task that was used to train the network – reconstruction / auto-encoding – against three baselines (trained in both single-category and multi-category variants): 1 3D-PointCapsNet [64], an auto-encoder for 3D point clouds that utilize a capsule architecture; 2 AtlasNetV2 [12], a state-of-the-art auto-encoder which utilizes a multi-head patch-based decoder; 3 AtlasNetV2 [12] with a spatial-transformer network (STN) aligner from PointNet [36], a baseline with canonicalization. We do not compare against [46], as unfortunately source code is not publicly available.
Quantitative analysis – Table 1. We achieve state-of-the-art performance in both the aligned and unaligned settings. The wider margin in the unaligned setup indicates tackling this more challenging scenario damages the performance of AtlasNetV2 [12] and 3D-PointCapsNet [64] much more than our method3. We also include a variant of AtlasNetV2 for which a STN (Spatial Transformer Network) is used to pre-align the point clouds [36], demonstrating how the simplest form of prealigner/canonicalizer is not sufficient.
Qualitative analysis – Figure 2. We illustrate our decomposition-based reconstruction of 3D point clouds, as well as the reconstructions of 3D-PointCapsNet [64] and AtlasNetV2 [12]. As shown, even in the unaligned setup, our method is able to provide semantically consistent capsule decompositions – e.g. the wings of the airplane have consistent colours, and when aligned in the canonical frame, the different airplane instances are well-aligned. Compared to AtlasNetV2 [12] and 3D-PointCapsNet [64], the reconstruction quality is also visibly improved: we better preserve details along the engines of the airplane, or the thin structures of the bench; note also that the decompositions
3Results in this table differ slightly from what is reported in the original papers as we use 1024 points to speed-up experiments throughout our paper. However, in the supplementary material the same trend holds regardless of the number of points, and match with what is reported in the original papers with 2500 points.
are semantically consistent in our examples. Results are better appreciated in our supplementary material, where we visualize the performance as we continuously traverse SE(3).

heading: 4.3 Canonicalization – Table 2

text:
We compare against three baselines: 1 Deep Closest Points [56], a deep learning-based pairwise point cloud registration method; 2 DeepGMR [63] a state-of-the-art pairwise registration method that decomposes clouds into Gaussian mixtures and utilizes Rigorously Rotation-Invariant (RRI) features [4]; 3 Compass [45] a concurrent work on learnt alignment/canonicalization. For all compared methods we use the official implementation. For DeepGMR we use both RRI and the typical XYZ coordinates as input. We also try our method with the RRI features, following DeepGMR’s training protocol and train for 100 epochs.
Metrics. To evaluate the canonicalization performance, we look into the stability of the canonicalization – the shakiness shown in the videos in our supplementary material– represented as the mean standard deviation of the rotations (mStd):
mStd = 1
n Xn i=1
r 1
m Xm j=1 (\(Rij ,Rimean))2 , (12)
where \ is the angular distance between two rotation matrices [18, 62, 63], Rij is the rotation matrix of the j-th instance of the i-th object in canonical frame, and Rimean is the mean rotation [21] of the i-th object. Note that with mStd we measure the stability of canonicalization with respect to rotations to accommodate for methods that do not deal with translation [45]. To allow for comparisons with pairwise registration methods, we also measure performance in terms of the RMSE metric [63, 66].
Quantitative analysis – Table 2. Compared to Compass [45], our method provides improved stability in canonicalization. This also provides an advantage in pairwise registration, delivering state-of-the-art results when XYZ-coordinates are used. Note that while pairwise methods can align
two sets of given point clouds, creating a canonical frame that simultaneously registers all point clouds is a non-trivial extension to the problem.
When RRI is used as input, our method is on par with DeepGMR [63], up to a level where registration is near perfect – alignment differences when errors are in the 10 4 ballpark are indiscernible. We note that the performance of Deep Closest Points [56] is not as good as reported in the original paper, as we uniformly draw rotations from SO(3). When a sub-portion of SO(3) is used, e.g. a quarter of what we are using, DCP performs relatively well (0.008 in the multi-class experiment). While curriculum learning could be used to enhance the performance of DCP, our technique does not need to rely on these more complex training techniques.
We further note that, while RRI delivers good registration performance, using RRI features cause the learnt canonicalization to fail – training becomes unstable. This hints that RRI features may be throwing away too much information to achieve transformation invariance. Our method using raw XYZ coordinates as input, on the other hand, provides comparable registration performance, and is able to do significantly more than just registration (i.e. classification, reconstruction).

heading: 4.4 Unsupervised classification – Table 3

text:
Beyond reconstruction and canonicalization, we evaluate the usefulness of our method via a classification task that is not related in any way to the losses used for training. We compute the features from the auto-encoding methods from Section 4.2 against those from our method (where we build features by combining pose with descriptors) to perform 13-way classification with two different techniques:
Table 3: Classification – Top-1 accuracy(%)
Method SVM K-Means
A lig
ne d 3D-PointCapsNet [64] 93.81 65.87
AtlasNetV2 [12] 94.07 61.66 Our method 94.21 69.82
U na
lig ne d 3D-PointCapsNet [64] 71.13 14.59 AtlasNetV2 [12] 64.85 17.12 AtlasNetV2 [12] w/ STN 78.55 20.03 Our method 87.33 43.04
• We train a supervised linear Support Vector Machine (SVM) on the extracted features [1, Ch. 7];
• We perform unsupervised K-Means clustering [1, Ch. 9] and then label each cluster via bipartite matching with the actual labels through the Hungarian algorithm.
Note the former provides an upper bound for unsupervised classification, while better performance on the latter implies that the learnt features are able to separate the classes into clusters that are compact (in an Euclidean sense).
Analysis of results – SVM. Note how our method provides best results in all cases, and when the dataset is not unaligned the difference is significant. This shows that, while 3D-PointCapsNet and AtlasNetV2 (with and without STN) are able to somewhat auto-encode point clouds in the unaligned setup, what they learn does not translate well to classification. However, the features learned with Canonical Capsules are more related to the semantics of the object, which helps classification.
Analysis of results – K-Means. The performance gap becomes wider when K-Means is used – even in the aligned case. This could mean that the features extracted by Canonical Capsules are better suited for other unsupervised tasks, having a feature space that is close to being Euclidean in terms of semantics. The difference is striking in the unaligned setup. We argue that these results emphasize the importance of the capsule framework – jointly learning the invariances and equivariances in the data – is cardinal to unsupervised learning [25, 24].

heading: 4.5 Ablation study

text:
To make the computational cost manageable, we perform all ablations with the airplane category (the category with most instances), and in the unaligned setup (unless otherwise noted). Please also see supplementary material for more ablation studies.
Losses – Table 4. We analyze the importance of each loss term, with the exception of Lrecon which is necessary for training. All losses beneficially contribute to reconstruction performance, but note how Lequiv, Llocalization and Lequilibrium affect it to a larger degree. By considering our canonicalization metric, we can motivate this outcome by observing that the method fails to perform canonicalization when these losses are not employed (i.e. training collapses).
Encoder architecture – Table 5. Our method can be used with any backbone, as our main contribution lies in the self-supervised canonicalization architecture. For completeness, we explore variants of our method using different back-bones: PointNet [36], PointNet++ [37], DGCNN [57], and ACNe [47]. Among these, the ACNe backbone performs best; note that all the variants in Table 5, regardless of backbone choice, significantly outperforms all other methods reported in Table 1.

heading: 5 Conclusions

text:
In this paper, we provide a self-supervised framework to train primary capsule decompositions for 3D point clouds. We rely on a Siamese setup that allows self-supervision and auto-encoding in canonical space, circumventing the customary need to train on pre-aligned datasets. Despite being trained in a self-supervised fashion, our representation achieves state-of-the-art performance across auto-encoding, canonicalization and classification tasks. These results are made possible by allowing the network to learn a canonical frame of reference. We interpret this result as giving our neural networks a mechanism to construct a “mental picture” of a given 3D object – so that downstream tasks are executed within an object-centric coordinate frame.
Future work. As many objects have natural symmetries that we do not consider at the moment [12], providing our canonicalizer a way to encode such a prior is likely to further improve the representation. We perform decomposition at a single level, and it would be interesting to investigate how to effectively engineer multi-level decompositions [51]; one way could be to over-decompose the input in a redundant fashion (with large K), and use a downstream layers that “selects” the decomposition heads to be used [6]. We would also like to extend our results to more “in-the-wild” 3D computer vision and understand whether learning object-centric representations is possible when incomplete (i.e., single view [35] or occluded) data is given in input, when an entire scene with potentially multiple objects is given [43], or where our measurement of the 3D world is a single 2D image [48], or by exploiting the persistence of objects in video [40].
Broader impact. While our work is exclusively on 3D shapes, and thus not immediately subject to any societal concerns, it enhances how Artificial Intelligence (AI) can understand and model 3D geometry. Thus, similar to how image recognition could be misused, one should be careful when extending the use of our method. In addition, while not subject to how the data itself is aligned, the learnt canonical frame of our method is still data-driven, thus subject to any data collection biases that may exist – canonicalization will favour shapes that appear more often. This should also be taken into account with care to prevent any biased decisions when utilizing our method within a decision making AI platform.

heading: Acknowledgements

text: This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant, NSERC Collaborative Research and Development Grant, Google, Compute Canada, and Advanced Research Computing at the University of British Columbia.

title: Canonical Capsules: Self-Supervised Capsules in Canonical Pose


summaries: This paper introduces a new method to construct conformal prediction intervals based on histograms of the conditional distribution of an outcome variable. Given a histogram of the conditional distribution of the outcome variable the method finds the shortest interval whose associated mass is no less than the desired coverage probability of the prediction interval. This yields intervals that automatically adapt to the skewness of the data. The intervals have provably correct marginal coverage in finite samples and correct conditional coverage and optimal length in large samples. A simulation study and numerical experiments on several benchmark data sets corroborate the theoretical results and demonstrate the advantage over several competing methods.

summaries: This paper proposes a new nonconformity-score based on conditional density estimates P(Y|X). The main idea is to approximate the optimal prediction interval based on the conditional density estimate. The proposed non-conformity score is experimentally assessed using the conformalization technique of split conformal. The overall method provably achieves marginal coverage, performs comparably to other methods for conditional coverage, and achieves the shortest prediction intervals on all datasets.

summaries: This paper proposes an extension to conformal prediction that adapts to skewed data, and can achieve better conditional coverage (and provably achieves conditional coverage asymptotically). Conformal prediction in general is a methodology for constructing confidence sets that output likely response candidates $\widehat{\mathcal{C}}(X) \subseteq \mathcal{Y}$ for and input $X$, rather than a single value. The goal is to ensure that $\widehat{\mathcal{C}}$ covers the true response variable, $Y$, with specifiably high probability. While many prior methods provably control *marginal* coverage, *conditional* coverage is a much harder (and more practically important) goal (albeit impossible in finite samples in the general case). This paper attacks this problem by developing a novel conformalization strategy that leverages calibrated estimates of the conditional density $Y \mid X$ to obtain approximate conditional coverage in finite samples (empirically), and asymptotic conditional coverage (theoretically).

summaries: This paper introduces a conformity score that aims to decrease average interval lengths and improve conditional coverage from a black-box estimate  $\hat{P}(y \mid x)$ for split conformal prediction. The method, conformal histogram regression (CHR), involves first binning the space of $Y$, resulting in a conditional histogram from which approximate oracle intervals can be computed. A nested sequence of these intervals is then created for a sequence of predictive miscoverage values $\tau$, where $\tau$ will be close to $\alpha$ if $\hat{P}(y \mid x)$ is a good estimate. The value of $\tau$ is selected through a conformity score, and the authors show their method obtains finite marginal coverage and asymptotic conditional coverage. The authors then demonstrate the method on a few examples.

summaries: The paper shows how to obtain *short* prediction intervals that approximately have the right conditional coverage. In order to do so, the method estimates the conditional distribution of a new label given its features, and then uses a conformal score based on the quantiles implied by such estimate to create the prediction intervals. The experiments show that the method outperforms other quantile-based conformal methods in terms of width while still approximately controlling conditional coverage. Some theoretical results that prove converge to the oracle (i.e., the shortest prediction interval) are also shown.

abstractText: This paper develops a conformal method to compute prediction intervals for nonparametric regression that can automatically adapt to skewed data. Leveraging black-box machine learning algorithms to estimate the conditional distribution of the outcome using histograms, it translates their output into the shortest prediction intervals with approximate conditional coverage. The resulting prediction intervals provably have marginal coverage in finite samples, while asymptotically achieving conditional coverage and optimal length if the black-box model is consistent. Numerical experiments with simulated and real data demonstrate improved performance compared to state-of-the-art alternatives, including conformalized quantile regression and other distributional conformal prediction approaches.

authors:
R.F. Barber
E. Candès
A. Ramdas
R.J. Tibshirani

title: Predictive inference with the jackknife+

venue: Annals of Statistics,

authors:
S. Bates
E. Candès
L. Lei
Y. Romano
M. Sesia

title: Testing for outliers with conformal p-values

authors:
M. Cauchois
S. Gupta
J.C. Duchi

title: Knowing what you know: valid and validated confidence sets in multiclass and multilabel prediction

venue: Journal of Machine Learning Research,

authors:
V. Chernozhukov
K. Wüthrich
Y. Zhu

title: Distributional conformal prediction

authors:
H.A. Chipman
E.I. George
R.E. McCulloch

title: BART: Bayesian additive regression trees

venue: Annals of Applied Statistics,

authors:
A. Clauset
C.R. Shalizi
M.E. Newman

title: Power-law distributions in empirical data

venue: SIAM review,

authors:
J.W. Cohen
S.B. Cohen
J.S. Banthin

title: The medical expenditure panel survey: a national information resource to support healthcare cost research and inform policy and practice

venue: Medical care, pages S44–S50,

authors:
N. Dalmasso
T. Pospisil
A.B. Lee
R. Izbicki
P.E. Freeman
A.I. Malz

title: Conditional density estimation tools in python and r with applications to photometric redshifts and likelihoodfree cosmological inference

venue: Astronomy and Computing,

authors:
R. Foygel Barber
E. Candès
A. Ramdas
R.J. Tibshirani

title: The limits of distribution-free conditional predictive inference

venue: Information and Inference: A Journal of the IMA,

authors:
C. Gupta
A.K. Kuchibhotla
A.K. Ramdas

title: Nested conformal prediction and quantile out-of-bag ensemble methods

title: Making learning more transparent using conformalized performance prediction

authors:
R. Izbicki
A.B. Lee

title: Nonparametric conditional density estimation in a high-dimensional regression setting

venue: Journal of Computational and Graphical Statistics,

authors:
R. Izbicki
G. Shimizu
R.B. Stern

title: CD-split: efficient conformal regions in high dimensions

authors:
R. Izbicki
G.T. Shimizu
R.B. Stern

title: Flexible distribution-free conditional predictive bands using density estimators

authors:
B. Kim
C. Xu
R.F. Barber

title: Predictive inference is free with the jackknife+-after-bootstrap

authors:
D. Kivaranovic
K.D. Johnson
H. Leeb

title: Adaptive, distribution-free prediction intervals for deep networks

venue: In International Conference on Artificial Intelligence and Statistics,

authors:
J. Lei
M. G’Sell
A. Rinaldo
R.J. Tibshirani
L. Wasserman

title: Distribution-free predictive inference for regression

venue: Journal of the American Statistical Association,

authors:
J. Lei
L. Wasserman

title: Distribution-free prediction bands for non-parametric regression

venue: Journal of the Royal Statistical Society: Series B (Statistical Methodology),

authors:
M. Magdon-Ismail
A. Atiya

title: Neural networks for density estimation

venue: In Proceedings of the 11th International Conference on Neural Information Processing Systems,

title: Quantile regression forests

authors:
S.J. Moon
J.-J. Jeon
J.S.H. Lee
Y. Kim

title: Learning multiple quantiles with neural networks

authors:
Y. Romano
R.F. Barber
C. Sabatti
E. Candès

title: With malice toward none: Assessing uncertainty via equalized coverage

venue: Harvard Data Science Review,

authors:
Y. Romano
E. Patterson
E. Candès

title: Conformalized quantile regression

authors:
Y. Romano
M. Sesia
E. Candès

title: Classification with valid and adaptive coverage

authors:
M. Sesia
E. Candès

title: A comparison of some conformal quantile regression methods. Stat

title: A quantile regression neural network approach to estimating the conditional density of multiperiod returns

venue: Journal of Forecasting,

title: Conditional validity of inductive conformal predictors

venue: In Asian conference on machine learning,

authors:
V. Vovk
A. Gammerman
G. Shafer

title: Algorithmic learning in a random world

authors:
Y. Yang
A.K. Kuchibhotla

title: Finite-sample efficient conformal prediction

venue: preprint at arXiv:2104.13871,

heading: 1 Introduction

heading: 1.1 Problem statement and motivation

text:
We consider the problem of predicting with confidence a response variable Y ∈ R given p features X ∈ Rp for a test point n+1, utilizing n pairs of observations {(X(i), Y (i))}ni=1 drawn exchangeably (e.g., i.i.d.) from some unknown distribution, and leveraging any machine-learning algorithm. Precisely, ∀α ∈ (0, 1), we seek a prediction interval Ĉn,α(Xn+1) ⊂ R for Yn+1 satisfying the following three criteria. First, Ĉn,α should have finite-sample marginal coverage at level 1− α,
P [ Yn+1 ∈ Ĉn,α(Xn+1) ] ≥ 1− α. (1)
Second, Ĉn,α should approximately have conditional coverage at level 1− α, P [ Yn+1 ∈ Ĉn,α(x) | Xn+1 = x ] ≥ 1− α, ∀x ∈ Rp, (2)
meaning it should approximate this objective in practice, and ideally achieve it asymptotically under suitable conditions in the limit of large sample sizes. Third, Ĉn,α should be as narrow as possible.
We tackle this challenge with conformal inference [24, 35], which allows one to convert the output of any black-box machine learning algorithm into prediction intervals with provable marginal coverage (1). The key idea of this framework is to compute a conformity score for each observation, measuring the discrepancy, according to some metric, between the true value of Y and that predicted
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
by the black-box model. The model fitted on the training data is then applied to hold-out calibration samples, producing a collection of conformity scores. As all data points are exchangeable, the empirical distribution of the calibration scores can be leveraged to make predictive inferences about the conformity score of a new test point. Finally, inverting the function defining the conformity scores yields a prediction set for the test Y . This framework can accommodate almost any choice of conformity scores, and in fact many different implementations have already been proposed to address our problem. However, it remains unclear how to implement a concrete method from this broad family that can lead to the most informative possible prediction intervals. Our contribution here is to develop a practical solution, following the three criteria defined above, that performs better compared to existing alternatives and is asymptotically optimal under certain assumptions.
It is worth emphasizing that constructing a short prediction interval with guaranteed coverage is a reasonable approach to quantify and communicate predictive uncertainty in regression problems, although it is of course not the only one. To name an alternative, one could compute a non-convex prediction set with analogous coverage [20], which might be more appropriate in some situations, but is also more easily confusing. For example, it could be informative for a physician to know that the future blood pressure of a patient with certain characteristics is predicted to be within the range [120,129] mmHg. However, it would not be more helpful to report instead the following non-convex region: [120, 120.012]∪ [120.015, 120.05]∪ [121, 122.7]∪ [123.1, 127.2]∪ [127.8, 129] mmHg. Indeed, in the second case it would not be clear (a) whether the multi-modal nature of that prediction is significant or a spurious consequence of overfitting, and (b) how the physician would act upon that prediction any differently than if it had been [120,129] mmHg. Therefore, we focus on prediction intervals in this paper because they are generally easier to interpret than arbitrary regions, and they are also less likely to convey a false sense of confidence.

heading: 1.2 Preview of conformal histogram regression

text:
Imagine an oracle with access to PY |X , the distribution of Y conditional on X , which leverages such information to construct optimal prediction intervals as follows. For simplicity, suppose PY |X has a continuous density f(y | x) with respect to the Lebesgue measure, although this could be relaxed with more involved notation. Then, the oracle interval for Yn+1 | Xn+1 = x would be:
Coracleα (x) = [ loracle1−α (x), u oracle 1−α (x) ] , (3)
where, for any τ ∈ (0, 1], loracleτ (x) and uoracleτ (x) are defined as:
[loracleτ (x), u oracle τ (x)] := argmin
(l,u)∈R2 : l≤u
{ u− l : ∫ u l f(y | x)dy ≥ τ } . (4)
This is the shortest interval with conditional coverage (2). If the solution to (4) is not unique (e.g., if f(· | x) is piece-wise constant), the oracle picks any solution at random. Of course, this is not a practical method because f is unknown. Therefore, we will approximate (4) by fitting a black-box model on the training data, and then use conformal prediction to construct an interval accounting for any possible estimation errors. Specifically, we replace f in (4) with a histogram approximation, hence why we call our method conformal histogram regression, or CHR. The output interval is then
Ĉn,α(x) = [ l̂τ̂ (x), ûτ̂ (x) ] , (5)
where l̂τ̂ (x) and ûτ̂ (x) approximate the analogous oracle quantities in (4). The value of τ̂ in (5) will be determined by suitable conformity scores evaluated on the hold-out data, and it may be larger than 1 − α if the model for f is not very accurate. However, if the fitted histogram is close to the true PY |X , the interval in (5) will resemble that of the oracle (3).
Figure 1 previews an application to toy data, comparing CHR to conformalized quantile regression (CQR) [30]; see Section 4.2 for more details. CHR finds the shortest interval such the corresponding area under the histogram is above τ , for any τ ∈ (0, 1], and then calibrates τ to guarantee marginal coverage above 1−α; this extracts more information from the model compared to CQR. For example, CHR adapts automatically to the skewness of Y | X , returning intervals delimited by the 0%–90% quantiles in this example, which are shorter than the symmetric ones (5%–95%) sought by CQR.

heading: 1.3 Related work

text:
This work is inspired by the conformity scores introduced by [31] for multi-class classification, the underlying idea of which can be repurposed here. Nonetheless, the extension to our problem involves several innovations. This paper connects [31] to other conformal methods for continuous responses [24, 25, 35], which sought objectives similar to ours by leveraging quantile regression [17, 23, 30, 32, 36] or non-parametric density estimation [10, 21], sometimes considering multi-modal prediction sets instead of intervals [20]. Our approach also exploits black-box models for distributional estimation; however, we introduce more efficient conformity scores.
We seek the shortest intervals with marginal coverage while approximating as well as possible conditional coverage, although the latter is impossible to guarantee in finite samples [16, 34]. The performances of prior approaches have been measured in terms of these criteria, yet others have not sought them as directly. Indeed, if the black-box model is consistent for PY |X , our method becomes asymptotically equivalent to the oracle (3)–(4), under some technical assumptions. This property does not hold for other existing methods because they tend to produce symmetric intervals, with fixed lower and upper miscoverage rates (the probabilities of the outcome being either below or above the output interval, respectively), which may be sub-optimal if the data have unknown skewness.

heading: 2 Methods

text: The proposed method consists of four main components: the estimation and binning of a conditional model for the outcome, the construction of a nested sequence of approximate oracle intervals based on the above, the computation of suitable conformity scores, and their conformal calibration.

heading: 2.1 Estimating conditional histograms

text:
We partition the domain of Y into m bins [bj−1, bj), for some sequence b0 < . . . < bm. With little loss of generality, assume Y is bounded: −C = b0 < Y < bm = C, for some C > 0. Then, we solve a discrete version of the problem stated in the introduction: we seek the smallest possible contiguous subset of bins with 1− α predictive coverage. If m is large and the bins are narrow, this problem is not very different from the original one, although it is more amenable to solution.
For simplicity, we present our method from a split-conformal perspective [24, 35]; extensions to other hold-out approaches [7, 22, 35] will be intuitive. Let Dtrain,Dcal ⊂ {1, . . . , n} denote any partition of the data into training and calibration subsets, respectively. Dtrain is used to train a black-box model for the conditional probabilities that Y is within any of the above bins: ∀j ∈ {1, . . . ,m},
πj(x) := P[Y ∈ [bj−1, bj) | X = x]. (6)
There exist many tools to approximate PY |X and obtain estimates π̂j(x) of πj(x), including quantile regression [27, 28, 33], Bayesian additive regression trees [11], or any other non-parametric conditional density estimator [14, 19, 26]. Our method can directly be applied with any of these models, but we found multiple quantile regression to work particularly well [17, 23, 30, 32], and therefore we will focus on it in this paper. Referring to Supplementary Section S1.1 for implementation details and information about the computational cost of the learning algorithm (which is comparable to that required by CQR [30]), we thus take these black-box estimates π̂j(x) as fixed henceforth.
Note that estimating conditional distributions is more challenging if the number of variables is larger. However, this is a fundamental difficulty of high-dimensional regression, not a particular limitation of the proposed CHR. Although our method utilizes conditional histograms learnt from the data, its performance is not directly measured in terms of how closely these resemble the true PY |X . Instead, as we shall see, CHR only needs to detect the possible skewness of Y | X and estimate reasonably well some lower and upper quantiles of this conditional distribution. Therefore, its estimation task is not much more difficult than that of CQR [30], as skewness is relatively easy to detect.

heading: 2.2 Constructing a nested sequence of approximate oracle intervals

text:
For any partition B = (b0, . . . , bm) of the domain of Y , let π = (π1, . . . , πm) be a unit-sum sequence, depending on x ∈ Rp; this may be seen as a histogram approximation of PY |X (6). For simplicity, assume all histogram bins have equal width, although this is unnecessary. Then, define the following bi-valued function S taking as input x ∈ Rp, π, τ ∈ (0, 1], and two intervals S−, S+ ⊆ {1, . . . ,m}:
S(x, π, S−, S+, τ) := argmin (l,u)∈{1,...,m}2 : l≤u |u− l| : u∑
j=l
πj(x) ≥ τ, S− ⊆ [l, u] ⊆ S+  . (7)
Above, it is implicitly understood we choose the value of (l, u) minimizing ∑u
j=l πj(x) among the feasible ones with minimal |u− l|, if the optimal solution would not otherwise be unique. Therefore, we can assume without loss of generality the solution to (7) is unique; if that is not the case, we can break the ties at random by adding a little noise to π. The problem in (7) can be solved at computational cost linear in the number of bins, and it is equivalent to the standard programming challenge of finding the smallest positive subarray whose sum is above a given threshold. Note that we will sometimes refer to intervals on the grid determined by B as either contiguous subsets of {1, . . . ,m} (e.g., S−) or as pairs of lower and upper endpoints (e.g., [l, u]). If S− = ∅ and S+ = {1, . . . ,m}, the expression in (7) computes the shortest possible interval with total mass above τ according to π(x). Further, if πj is the mass in the j-th bin according to the true PY |X , then S(x, π, ∅, {1, . . . ,m}, 1 − α) is the discretized version of the oracle interval (3)– (4). In general, the optimization in (7) involves the additional nesting constraint that the output S must satisfy S− ⊆ S ⊆ S+, which will be needed to guarantee our method has valid marginal coverage [17]. Intuitively, it is helpful to work with a nested sequence because this ensures the prediction intervals are monotone increasing in τ , essentially reducing the calibration problem to that of selecting the appropriate value of τ that yields the desired marginal coverage. Note that the inequality in (7) involving τ may not be binding at the optimal solution due to the discrete nature of the optimization problem. However, the above oracle can be easily modified by introducing some suitable randomization in order to obtain valid prediction intervals that are even tighter on average, as explained in Supplementary Section S1.2.
As π̂ may be an inaccurate estimate of PY |X , we cannot simply plug it into the oracle in (7) and expect valid coverage. However, for any approximate conditional histogram π̂, we can define a nested sequence [17] of (randomized) sub-intervals of B, for different values of τ ranging from 0 to 1. Then, we calibrate τ to obtain the desired 1− α marginal coverage. Precisely, consider an increasing scalar sequence τt = t/T , for t ∈ {0, . . . , T} with some T ∈ N, and define a corresponding growing sequence of subsets St ⊆ {1, . . . ,m} as follows. First, fix any starting point t̄ ∈ {0, . . . , T} and define St̄ by applying (7) without the nesting constraints (with S− = ∅ and S+ = {1, . . . ,m}):
St̄ := S(x, π, ∅, {1, . . . ,m}, τt̄), (8) Note the explicit dependence on x and π of the left-hand-side above is omitted for simplicity, although it is important to keep in mind that St̄ does of course depend on these quantities. Figure 2 (second row) visualizes the construction of St̄ in a toy example with τt̄ = 0.9.
Having computed the initial interval St for t = t̄, we recursively extend the definition to the wider intervals indexed by t = t̄+ 1, . . . , T as follows:
St := S(x, π, St−1, {1, . . . ,m}, τt). (9) See the top row of Figure 2 for a schematic of this step. Similarly, the narrower intervals St indexed by t = t̄− 1, t̄− 2, . . . 0 are defined recursively as:
St := S(x, π, ∅, St+1, τt). (10) See the bottom row of Figure 2 for a schematic of this step. As a result of this construction, the sequence of intervals {St}Tt=0 is nested regardless of the starting point t̄ in (8), as previewed in Figure 2. However, different choices of t̄ may lead to different sequences, any of which allows us to obtain provable marginal coverage, as discussed next. As our goal is to approximate the oracle in (3)–(4) accurately, the most intuitive choice is to pick t̄ such that τt̄ ≈ 1 − α. A more involved randomized version of this construction, inspired by the more powerful randomized oracle, is discussed in Supplementary Section S1.2. Note that the randomized version of the nested prediction intervals will be the one applied throughout this paper and, with a slight overload of notation, we will simply refer to it as {St}Tt=1. Note also that we will highlight the dependence of this sequence on x and π by writing it as St(x, π). Further, as we work henceforth with the randomized versions of these prediction intervals (described in Supplementary Section S1.2), we will refer to them as St(x, ε, π), where ε is a uniform random variable in [0, 1], independent of everything else.

heading: 2.3 Computing conformity scores and calibrating prediction intervals

text:
Given any sequence of nested sets St(x, ε, π), we define the following conformity score function E: E(x, y, ε, π) := min {t ∈ {0, . . . , T} : y ∈ St(x, ε, π)} . (11)
In words, this computes the smallest index t such that St(x, ε, π) contains y, as in [17, 31]. Equivalently, one can think of these scores as indicating the smallest value of the nominal coverage τ in (7) necessary to ensure the observed Y is contained in the prediction interval. Our method evaluates (11) on all calibration samples (Xi, Yi) using the π̂ learnt on the training data; for each i ∈ Dcal, we generate εi ∼ Unif(0, 1) and store
Ei = E(Xi, Yi, εi, π̂).
Then, we compute prediction intervals for Yn+1 by looking at the nested sequence in (8)–(10) corresponding to the new Xn+1 and selecting the interval indexed by the 1− α quantile (roughly) of {Ei}i∈Dcal . The procedure is outlined in Algorithm 1. Note that the only computationally expensive component of this method is the estimation of the conditional histograms (see Supplementary Section S1.1 for details); the construction of the nested prediction intervals and the evaluation of the conformity scores have negligible cost because the optimization problem in (7) is an easy one.
It may be helpful to point out that, if π provides an accurate representation of the true conditional distribution of Y | X , then the above conformity scores are uniformly distributed [31]. In that ideal case, no calibration is needed and indeed our method simply reduces to applying (7) with τ = 0.9 to construct prediction intervals with 90% coverage. In practice, however, π can only be a possibly inaccurate estimate of PY |X (hence why we will refer to it as π̂ from now on), which means that the distribution of the conformity scores may not be uniform and the conformal calibration is necessary to obtain valid coverage.
The next result states that the output of our method has valid marginal coverage, regardless of the accuracy of π̂. The proof relies on the sequence St being nested; from there, coverage follows from the results of [17, 31]; see Supplementary Section S1.1.
Algorithm 1: CHR with split-conformal calibration 1 Input: data {(Xi, Yi)}ni=1, Xn+1, partition B of the domain of Y into m equal-sized bins, level
α ∈ (0, 1), resolution T for the conformity scores, starting index t̄ for recursive definition of conformity scores, black-box algorithm for estimating conditional distributions.
2 Randomly split the training data into two subsets, Dtrain,Dcal. 3 Sample εi ∼ Uniform(0, 1) for all i ∈ {1, . . . , n+ 1}, independently of everything else. 4 Using the data in Dtrain, train any estimate π̂ of the mass of Y | X for each bin in B (6); see
Supplementary Section S1.1 for a concrete approach based on quantile regression. 5 Compute Ei = E(Xi, Yi, εi, π̂) for each i ∈ Dcal, with the function E defined in (11). 6 Compute t̂ = Q̂1−α({Ei}i∈Dcal) as the ⌈(1− α)(1 + |Dcal|)⌉th smallest value in {Ei}i∈Dcal . 7 Select the t̂-th element from {St(Xn+1, εn+1, π̂)}Tt=0, defined in (8)–(10):
Ĉscn,α(Xn+1) = St̂(Xn+1, εn+1, π̂).
8 Output: A prediction interval Ĉscn,α(Xn+1) for Yn+1.
Theorem 1 (Marginal coverage). If (Xi, Yi), for i ∈ {1, . . . , n + 1}, are exchangeable, then the output of Algorithm 1 satisfies:
P [ Yn+1 ∈ Ĉscn,α(Xn+1) ] ≥ 1− α. (12)
Note that Theorem 1 provides only a lower bound; a nearly matching upper bound on the marginal coverage can be generally established for split-conformal inference if the conformity scores are almost-surely distinct [24, 30, 35]. Although the CHR scores (11) are discrete, our experiments will show the coverage is tight as long as the resolution T is not too small.

heading: 3 Asymptotic analysis

text:
We prove here that the prediction intervals computed by CHR (Algorithm 1) are asymptotically equivalent, as n → ∞, to those of the oracle from (3)–(4), if the model π̂ is consistent for PY |X and a few other technical conditions are met. In particular, we analyze a slightly modified version of Algorithm 1 in which there is no randomization; this is theoretically more amenable and equivalent in spirit, although it may yield wider intervals in finite samples. Our theory relies on the additional Assumptions 1–5, explained below and stated formally in Supplementary Section S2.
1. The samples are i.i.d., which is stronger than exchangeability; this is the key to our concentration results.
2. The black-box model estimates PY |X consistently, in a sense analogous to that in [24, 32]. This assumption is crucial and may be practically difficult to validate in practice, but it
can be justified by existing consistency results available for some models under suitable conditions, such as random forests [27]. Further, the resolution m of the partition of the Y domain should grow with n at a certain rate, and the resolution T of the scores Ei in (11) should grow as Tn = n.
3. The true PY |X is continuous and with bounded density within a finite domain. This assumption is technical and could be relaxed with more work.
4. The true PY |X is unimodal. This assumption is also technical and could be relaxed with more work.
5. The estimated conditional histogram π̂ preserves the boundedness and unimodality of PY |X ; this assumption may be unnecessary but it is convenient and quite innocuous at this point given Assumptions 2–4.
For simplicity, we assume the number of observations is 2n, the test point is (X2n+1, Y2n+1), and Dtrain = Dcal = n, although different relative sample sizes would yield the same results. Theorem 2 (Asymptotic conditional coverage and optimality). ∀α ∈ (0, 1], let Ĉscn,α(X2n+1) denote the prediction interval for Y2n+1 computed by Algorithm 1 at level 1 − α without randomization. Under Assumptions 1–5, Ĉscn,α(X2n+1) is asymptotically equivalent, as n → ∞, to Coracleα (X2n+1), the output of the oracle (3)–(4). In particular, the following two properties hold.
(i) Asymptotic oracle length. For some sequences γn → 0 and ξn → 0 as n → ∞, P [ |Ĉscn,α(X2n+1)| ≤ |Coracleα (X2n+1)|+ γn ] ≥ 1− ξn.
(ii) Asymptotic conditional coverage. For some sequences ϵn → 0 and ζn → 0 as n → ∞, P [ P [ Y ∈ Ĉscn,α(X2n+1) | X2n+1 ] ≥ 1− α− ϵn ] ≥ 1− ζn.
Theorem 2 is similar to results in [24] and [32] about the efficiency of earlier approaches to conformal regression, including CQR [32]. However, the increased flexibility of our method is reflected by the oracle in Theorem 2, which is stronger than those in [24, 32]. In fact, the oracle in [24] does not have conditional coverage, and that in [32] produces wider prediction intervals with constant lower and upper miscoverage rates. Other conformal methods based on non-parametric density estimation [10, 21] are not as efficient as CHR, in the sense that Theorem 2 does not hold for them.

heading: 4 Numerical experiments

heading: 4.1 Software implementation

text: A Python implementation of CHR is available online at https://github.com/msesia/chr, along with code to reproduce the following numerical experiments. This software divides the domain of Y into a desired number of bins with equal sizes, depending on the range of values observed in the training data; we use 100 bins for the synthetic data and 1000 for the real data. Then, we estimate the conditional histograms π̂ using different black-box quantile regression models [27, 33], with a grid of quantiles ranging from 1% to 99%; see Supplementary Section S1.1. Our software also supports Bayesian additive regression trees [11] and could easily accommodate other alternatives. For simplicity, we apply CHR and other benchmark methods by assigning equal numbers of samples to the training and calibration sets; this ensures all comparisons are fair, although different options may lead to even shorter intervals [32]. See [8] for a rigorous discussion of how this choice affects the variability of the coverage conditional on the calibration data, which is an issue we do not explore here. See Supplementary Section S3 for details about how the models are trained, and information about the necessary computational resources.

heading: 4.2 Synthetic data

text:
We simulate a synthetic data set with a one-dimensional feature X and a continuous response Y , from the same distribution previewed in Figure 1, which is similar to that utilized in [30] to present CQR. Our method is applied to 2000 independent observations from this distribution, using the first 1000 of
them for training a deep quantile regression model, and the remaining ones for calibration. Figure 1 visualizes the resulting prediction bands for independent test data, comparing them to the analogous quantities output by CQR. Both methods are based on the same neural network and guarantee 90% marginal coverage, but ours leads to narrower intervals. Indeed, the advantage of CHR is that it can extract information from all conditional quantiles estimated by the base model and then automatically adapt to the estimated data distribution. By contrast, CQR [30] can only leverage a pre-specified lower and upper quantile (e.g., 5% and 95% in this example), and is therefore not adaptive to skewness.
Figure 3 (a) summarizes the performance of CHR over 100 experiments based on independent data sets, as a function of the sample size. We evaluate the marginal coverage, approximate the worst-slab conditional coverage [9] as in [31], and compute the average interval width. We consider two benchmarks in addition to CQR [30]: distributional conformal prediction (DCP) [10] and DistSplit [21]. To facilitate the comparisons, all methods have the same base model. (We also applied DistSplit as implemented by [21], with a different base model, but the version presented here performs better.) These results show CHR leads to the shortest prediction intervals, while simultaneously achieving the highest conditional coverage. Compatibly with Theorem 2, the output of CHR becomes roughly equivalent to that of the omniscient oracle as the sample size grows; the latter can be implemented exactly here because we know the true data generating process.
Supplementary Figure S1 compares the performance of CHR in these experiments to that of naive uncalibrated 90% prediction intervals based on the same deep neural network regression model, and obtained by simply plugging π̂ into the oracle in (7), with τ = 0.9. Unsurprisingly, the naive prediction intervals do not generally have the desired marginal coverage; in this case, they tend to be too narrow if the sample size is small and too wide if the sample size is large. Although the lack of coverage for the uncalibrated intervals is not very pronounced here because this black-box model is relatively accurate even with n = 100, such naive approach can yield arbitrarily low coverage in general, especially if the learning task is more difficult (e.g., for high-dimensional X), and is thus not reliable.
Figure 3 (b) shows analogous results from experiments in which we fix the sample size to 1000 and vary instead the skewness of the data distribution. Precisely, we flip a biased coin for each data point and transform Y into −Y if it lands heads, varying the coin bias as a control parameter. At one end of this spectrum, we recover the same skewed data distribution as in Figure 3 (a); at the other end, we have a symmetric PY |X . Our results are reported as a function of the expected skewness, defined
as E[(Y − µ(X))3/σ3(X)], where µ(X) and σ(X) are the mean and standard deviation of Y | X , respectively. These experiments show all methods are equivalent in terms of interval length if PY |X is symmetric (skewness close to 0), while CHR can be much more powerful if PY |X is skewed.

heading: 4.3 Real data

text:
We apply CHR to the following seven public-domain data sets also considered in [30]: physicochemical properties of protein tertiary structure (bio) [6], blog feedback (blog) [1], Facebook comment volume [2], variants one (fb1) and two (fb2), from the UCI Machine Learning Repository [15]; and medical expenditure panel survey number 19 (meps19) [3], number 20 (meps20) [4], and number 21 (meps21) [5], from [13]. We refer to [30] for more details about these data. As in the previous section, we would like to compare CHR to CQR, DistSplit, and DCP. However, as DCP [10] is unstable on all but one of these data sets, sometimes leading to very wide intervals, we replace it instead with a new hybrid benchmark that we call DCP-CQR. This improves the stability of DCP by combining it with CQR [30], as explained in Supplementary Section S1.3. This limitation of DCP may be explained by noting the method needs to learn a reasonably accurate approximation of the full conditional distribution of Y | X , and its performance is particularly sensitive to the estimation of the tails, which is most difficult; see Supplementary Section S1.3 for more details. By contrast, CHR is robust because it only needs to estimate a histogram with relatively few bins—a much easier statistical task—and then it specifically focuses on finding the shortest intervals containing high probability mass. We apply all methods, including our CHR, based on the same deep quantile regression model. Their performances are evaluated as in the previous section, averaging over 100 independent experiments per data set. In each experiment, 2000 samples are used for training, 2000 for calibration, and the remaining ones for training. All features are standardized to have zero mean and unit variance. The nominal coverage rate is 90%.
Figure 4 shows the distribution of the conditional coverage and interval width corresponding to different methods, separately for each data set. To simplify the plots by using a shared vertical axis, the widths of the prediction intervals are scaled, separately for each data set, so that the smallest one is always equal to one. Marginal coverage is omitted here because all methods provably control it; however, it can be found in Supplementary Table S1. All methods perform well in terms of worst-slab conditional coverage. CHR outperforms the others in terms of statistical efficiency, as it consistently leads to the shortest intervals. CQR and DistSplit are comparable to each other, while DCP-CQR sometimes outputs wider intervals. Supplementary Figure S2 shows that analogous results
are obtained if we utilized a random forest model instead of a deep neural network. Supplementary Table S1 summarizes these results in more detail, including marginal coverage and the omitted performance of the original DCP. Finally, Supplementary Figure S3 compares the performance of CHR in these experiments with real data to that of naive uncalibrated 90% prediction intervals based on the same deep neural network regression model, as in Figure S1. These results show that the naive prediction intervals do not generally have the desired marginal coverage; in some cases they are too narrow, and in others they are too wide.

heading: 5 Conclusions

text:
This paper developed CHR, a non-parametric regression method based on novel conformity scores leading to shorter prediction intervals with coverage, and enjoying stronger asymptotic efficiency, compared to the state-of-the-art alternatives. Of course, real data are finite and our theory relies on assumptions which may be difficult to validate; nonetheless, it is a sanity check and it provides an informative comparison. Further, the experiments confirm CHR performs well in practice.
The ability of CHR to automatically adapt to unknown skewness may prove useful in practice, as empirical data often follow distributions with power-law tails [12]. Indeed, the data sets analyzed in Section 4.3 tend to have highly skewed outcomes with many observations equal to zero. At the same time, a limitation of CHR is that it does not rigorously control the lower and upper miscoverage rates, which may be important for some applications; if that is the case, the modified version of CQR proposed by [30] would be a better choice. Note that the standard implementations of CQR and of the other benchmarks [10, 21] considered in this paper are not guaranteed to separately control the lower and upper miscoverage rates. In any case, users of our method could naturally obtain approximations of the lower and upper miscoverage rates for any prediction interval by looking at the underlying conditional histograms, although these estimates are of course not calibrated in finite samples.
Algorithm S1 in Supplementary Section S1.4 extends CHR to accommodate cross-validation+ [7], which is often more powerful, and computationally expensive, compared to the split-conformal approach presented in this paper. The strategy is the same as that followed by [31] in the classification setting, although it requires an extra step, in which the standard cross-validation+ prediction set [7] is replaced by its convex hull to guarantee the final output is an interval [17]. Supplementary Theorem S1 establishes that Algorithm S1 leads to marginal coverage above 1− 2α, applying the more general theory from [17]. Finally, a possible directions for future research may involve the extension of our method to deal with multi-dimensional responses Y .
From a broader perspective, this paper falls within a rapidly growing body of works focusing on improving the interpretability and statistical reliability of machine learning algorithms. Prediction intervals with marginal coverage provide a convenient way of communicating uncertainty about the accuracy of any machine learning black-box, which is important to increase their reliability, to ensure their fairness [29], and to facilitate their acceptance. Furthermore, conformal prediction intervals provide a principled metric by which to compare different machine learning algorithms [18].
Acknowledgments and Disclosure of Funding
The authors are grateful to Stephen Bates, Emmanuel Candès, and Wenguang Sun for providing insightful comments about an earlier version of this manuscript. M.S. thanks the center for Advanced Research Computing at the University of Southern California for providing computing resources. Y.R. was supported by the Israel Science Foundation (grant 729/21). Y.R. also thanks the Career Advancement Fellowship, Technion, for providing research support.

title: Conformal Prediction using Conditional Histograms