Technical Summary: Super-Level-Set Regression (SLS)

1. Problem Statement

In multivariate regression, accurately characterizing the conditional distribution of a response variable $Y \in \mathbb{R}^d$ given features $X$ is often computationally infeasible and statistically challenging. While full density estimation is possible, many decision-making tasks only require identifying the highest density regions (HDRs) that capture a specific probability mass $\tau$. These regions, known as conditional level sets, are traditionally extracted via a "plug-in" approach: first estimating the full conditional density $p(y|X)$ and then thresholding it.

This indirect method suffers from three primary limitations:

1. Inefficiency: It forces the model to expend capacity learning irrelevant parts of the distribution.
2. Sensitivity: It is highly sensitive to estimation errors in the density, which propagate to the thresholding step.
3. Rigidity: Direct estimation strategies often rely on restrictive prior assumptions regarding the shape, location, or topology of the data (e.g., convexity or unimodality).

Furthermore, formulating a direct optimization to minimize the volume of the prediction region subject to a conditional coverage constraint is difficult. The objective function (volume) is coupled with the conditional quantile of the model's own estimation error. Differentiating through this instance-dependent quantile is notoriously difficult, as existing techniques for quantile-constrained optimization (e.g., differentiable sorting, soft Lagrangian penalties) are designed for marginal or batch-level quantiles, not those varying dynamically with covariates $X$.

2. Methodology: Super-Level-Set Regression (SLS)

The authors introduce Super-Level-Set Regression (SLS), a framework that bypasses full density estimation to directly parameterize and optimize the geometric boundaries of conditional level sets.

Core Formulation

The goal is to find a mapping $\mathcal{A}(\cdot)$ that defines the conditional super-level set:
$$\mathcal{A}(X) = \{y \in \mathcal{Y} : p(y|X) \geq t_\tau(X)\}$$
where $t_\tau(X)$ is the threshold ensuring $P(Y \in \mathcal{A}(X) | X) \geq \tau$. Theoretically, the minimum-volume region satisfying this is the solution to:
$$\inf_{\mathcal{B}(X) \subset \mathcal{Y}} \text{Vol}(\mathcal{B}(X)) \quad \text{s.t.} \quad P(Y \in \mathcal{B}(X) | X) \geq \tau$$

Instead of estimating $p(y|X)$, SLS searches over a hypothesis space of real-valued functions $\mathcal{G}$ (frontier functions) to parameterize the region as a sub-level set:
$$\mathcal{C}_{G,q}(X) = \{y \in \mathcal{Y} : G(X, y) \leq q(X)\}$$
The optimization problem becomes minimizing the expected volume $\mathbb{E}_X[\text{Vol}_G(q(X), X)]$ subject to $q(X) \geq \text{Quantile}_\tau(P_{G(X,Y)}|X)$.

Solving the Implicit Coupling

The primary challenge is that the volume objective and the quantile constraint are implicitly coupled, making direct optimization intractable. To resolve this, the authors propose a surrogate objective that averages the volume functional over a shrinking probability neighborhood around the target quantile $\tau$.

Let $\varphi(n)$ and $\psi(n)$ be positive sequences converging to 0 as training steps $n \to \infty$. The surrogate objective minimizes:
$$J_n(G) = \mathbb{E}_{X,Y} \left[ K_n(X, G(X,Y)) \cdot \text{Vol}_G(G(X,Y), X) \right]$$
where $K_n$ is an indicator function restricting the expectation to the interval $[\tau - \varphi(n), \tau + \psi(n)]$.

Convergence: Proposition 3.1 establishes that as $n \to \infty$, this surrogate objective converges uniformly to the true constrained objective. This allows the use of an alternating optimization scheme:

1. Initialization: A "warm-up" phase where the window is wide (or the weighting is uniform), minimizing the average volume across all quantiles to establish a robust functional representation.
2. Alternating Updates:
   • Update the frontier function $G$ (parameters $\theta$) to minimize the volume within the current shrinking window.
   • Update the quantile estimator $q$ (parameters $\omega$) using the pinball loss to track the conditional quantiles of the current $G$.
3. Shrinking Window: The window boundaries $\varphi(n)$ and $\psi(n)$ decay over time, focusing the optimization strictly on the target $\tau$-quantile.

Frontier Functions

To handle complex geometries (multimodal, disjoint, non-convex), the paper proposes specific frontier architectures:

• Flow-based Mahalanobis Frontier: Uses a conditional volume-preserving normalizing flow $T_\phi$ combined with a learnable location $\mu_\theta(X)$ and shape matrix $L_\theta(X)$. The volume of the sub-level set admits a closed-form expression proportional to $q(X)^{d/2} / \det(L_\theta(X))$. This allows the model to capture asymmetric, star-shaped, or elliptical densities.
• Union of Flows: To model disjoint or highly multimodal regions, the method aggregates multiple flow components using a softmax-weighted average (approximating a hard minimum). This allows the prediction set to adapt to disconnected topologies without requiring extreme distortions in a single mapping.

3. Key Contributions

1. General Optimization Framework: The authors provide a theoretical framework (Proposition 3.1) capable of minimizing objective functions implicitly coupled with the conditional quantiles of the model's own predictions.
2. SLS Regression: A novel approach that focuses exclusively on learning specific conditional level sets, bypassing the intermediate step of full conditional density estimation and removing the need for restrictive prior assumptions about data shape.
3. End-to-End Geometric Optimization: A framework where the geometric boundary is directly parameterized and optimized, offering a new perspective on multivariate conditional quantile regression.
4. Adaptive Flexibility: By designing flexible parametric frontier functions (including unions of flows), the method adapts to the feature vector $X$, capturing complex structures such as disjoint and multimodal shapes.

4. Experimental Results

The authors evaluate SLS on synthetic and real-world datasets, comparing it against conditional density estimation (CDE) baselines and various conformal prediction strategies (e.g., C-HDR, PCP, L-CP).

• Synthetic Experiments:
  • SLS successfully recovers complex geometric shapes, including asymmetric, star-shaped densities and disjoint multimodal distributions.
  • In scenarios with outliers (e.g., Gaussian noise with 10% outliers), SLS learns to ignore the outliers to capture the highest density region, whereas standard maximum likelihood estimation (NLL) is heavily skewed by them.
  • The shrinking window mechanism is shown to be critical for correcting sub-optimal centering caused by the initialization phase in asymmetric distributions.
• Real-World Datasets:
  • Evaluated on six datasets (e.g., Bias, CASP, House, Taxi) with target coverages of 50% and 90%.
  • SLS performs competitively with CDE, operating directly at the Pareto frontier of the trade-off between region size (volume) and conditional coverage deviation.
  • While some baselines produce smaller sets, they often do so at the cost of significantly degraded conditional coverage. SLS maintains a better balance, achieving lower conditional coverage deviation for comparable volumes.

5. Significance and Claims

The paper claims that SLS represents a fundamental shift in perspective for multivariate conditional quantile regression. By replacing the restrictive assumptions of density-first methods with a direct geometric optimization strategy, SLS:

• Native Handling of Complexity: It natively captures complex, multimodal, and disjoint conditional structures without needing to model the full distribution.
• Direct Optimization: It resolves the implicit coupling between volume minimization and conditional quantile constraints through a novel surrogate objective, avoiding the need for differentiable sorting or marginal relaxations.
• Robustness: It is robust to outliers and misspecified distributional assumptions (e.g., assuming ellipticity when the true distribution is asymmetric) because it optimizes the region directly rather than fitting a global density.

The authors acknowledge that exact distribution-free guarantees for conditional coverage are theoretically impossible (Foygel Barber et al., 2021); thus, their method approximates this in practice, often supplemented by a post-hoc conformal calibration step to guarantee marginal validity. However, the primary contribution is the ability to learn highly adaptive, minimum-volume regions directly, offering a more efficient and flexible alternative to traditional plug-in density estimation.