Non-Asymptotic Uncertainty Quantification in High-Dimensional Learning

We sincerely thank you for your thorough and constructive feedback. Your insights have highlighted important areas for improvement in our paper's clarity and presentation. We would like to emphasize our novel contribution (see general rebuttal) and we kindly ask you to check the other reviews to increase your appreciation for the paper. Regarding the notation we plan to include a comprehensive "notation dictionary" in our paper, bridging the gap between statistics, signal processing, and machine learning communities. Furthermore, we appreciate your point about uncertainty quantification and will incorporate a detailed discussion on aleatoric versus epistemic uncertainty. We want to assure you that we have carefully addressed all of your concerns in our responses below, including clarifying acronyms, improving figure presentation, and elucidating the relationships between variables in our equations. That being said, we kindly ask you to reconsider your score and raise it. We respectfully would like to suggest that a difference in notation understanding or a few small points about the presentation (one figure which is poorly presented or the lack of some acronyms definition) may not be sufficient grounds for rejecting a paper. We are very committed to improving the presentation to make it accessible to a wider audience and will be happy to answer any further questions.

• W1. The paper uses non-standard notation, making it difficult to read. In Theorem 1, $b$ seems to represent the target, which is typically denoted as $y$. Additionally, the relationship between the matrix $A$ and the vectors $x$ is unclear. The l-1 norm in equation (1) is applied to $x$, but $x$ is also referred to as IID data in Theorem 1. Generally, the l-1 norm is used to penalize weights, commonly denoted by $w$, $\theta$, or $\beta$, rather than the input data for LASSO regression.

  • Thanks for pointing this out. We will make extensive comments in the final version to make the notation very clear for the different communities reading the paper. We will substitute $b$ with $y$ as it is more common in statistics and machine learning literature. Regarding the relationship between $A$ and $x$, let us describe the problem setting. The linear model: $y = Ax^0+\varepsilon$ consists of a matrix $A\in\mathbb{C}^{m\times N}$ and one ground truth vector $x^0\in\mathbb{C}^N$. Such notation is very common in the inverse problems literature. One of our goals is to recover $x^0$ when we know $y, A$ and the distribution of $\varepsilon$. If $x^0$ is sparse (e.g., radar images or angiographie images), one of the most common techniques in machine learning is the LASSO, i.e., we solve $\arg\min_x\Vert Ax - y\Vert_2^2 + \Vert x\Vert_1$. The $\ell_1$ norm is usually used to penalize nonzero entries, i.e., the regression coefficients, regardless of whether the vector represents NN weights or, as in our case, the ground truth vector, e.g., an image obtained from physical measurements. This regularization induces sparsity. If $x^0$ is not sparse, we use deep learning to obtain a reconstruction $x$ since some architectures, like the example we've given in the paper, are state-of-the-art for certain inverse problems. In this case, the training data for the deep learning model consists of i.i.d. feature vectors $x^{(1)},..., x^{(l)}$, following the same distribution as $x^0$ (e.g., could be similar magnetic resonance images like images from the same part of the body but from different patients) and the corresponding target vectors $y^{(i)} = A x^{(i)} +\varepsilon^{(i)}.$ We train the network to obtain a function $f$ that approximates $f(y^{(i)})\approx x^{(i)}$ and has the property $f(y^0)\approx x^0$. The point of our work is that for the first time as far as we can tell, we can quantify uncertainty componentwise very efficiently in a non-asymptotic way associated with the reconstruction of a given ground truth. The method not only comes with theoretical guarantees, but it is computationally cheap to implement since you don't need to recalculate the solution or retrain the model. We will be happy to provide further explanation here.

• W2. The paper does not clearly state the type of uncertainty being quantified, which could be clarified by addressing the first issue.

  • Thanks for this important question. We will add a discussion to the final version. We argue that we quantify the entire estimation uncertainty. One of the nice aspects of our method is that the decomposition of the estimation error into a Gaussian term and a remainder term allows for handling both types of uncertainty almost separately. With the Gaussian term $W$, we quantify the aleatoric uncertainty from the inherent measurement noise. The remainder term $R$ handles the epistemic uncertainty, which we quantify using a purely data-driven approach since for two different backward models (e.g., two different neural networks), it can be used to compare the estimation error of both models with respect to the ground truth. In this sense, our technique is rather a more general inferential uncertainty method.

• W3. Some acronyms are not defined (e.g., MR, ISTA, LASSO).

  • Thank you for pointing this out. The acronyms are all added to the paper. MR stands for Magnetic Resonance, the medical imaging modality, ISTA for terative Shrinkage Thresholding Algorithm, the most famous algorithm to solve the LASSO, and LASSO for Least Absolute Shrinkage and Selection Operator, the most important method for high-dimensional sparse regression.

• W4. Figure 1 is poorly presented. The images are very small with excessive white space in between, forcing readers to zoom in significantly. As a result, the caption becomes difficult to read.

  • Thanks, you are very correct. We will change the figures and labels/caption. Please check the pdf to see some of the new ones. Also, we will have one extra page which will allow us to have larger figures.

We are particularly grateful to the reviewer for raising the insightful point regarding conformal prediction and we will add some clarification about the difference between this technique and our work. While the two methods address aspects of total uncertainty, including both epistemic and aleatoric components, we carefully explain the fundamental differences between our approach and conformal prediction in our response below. We provide a detailed discussion on how our method complements rather than competes with conformal prediction, potentially opening avenues for future research that could bridge these approaches. We believe this clarification significantly enhances the positioning and contribution of our work within the broader landscape of uncertainty quantification in ML, and we would appreciate it if the reviewer could take this into account to raise our score.

• W1. The main weakness is that the proposed method is a competitor to conformal prediction. However, there is no comparison of these methods or even mention of this....

  • Thanks for the very relevant comment. We will add a discussion to the final version of the paper. Indeed, both methods produce a confidence/prediction interval for the output/prediction. However, they are inherently different approaches. Conformal prediction shines in generating prediction intervals for new observations (images) given the previous ones. On the other hand, the debiased LASSO produces precise confidence intervals for individual regression coefficients (pixels of a given image). The debiasing step corrects the bias made by the model; it is also particularly suited for problems in which the design matrix $A$, as well as the noise distribution, is known, e.g., inverse problems. Our method relies on additional samples (images $X_k$ and data $y_k$) only to calculate the distribution of the rest term $R$. This step of our method is always, as one can see in the experiments, the smallest portion of the error. Also, we don't need any calibration step that is computationally expensive in a regression setting. The distribution of the rest term is estimated by a random set of size $l$ chosen from the available data called the estimation data set. In this way, we can produce rigorous confidence intervals for each pixel of, for example, one single image.

    In contrast, conformal prediction solely relies on data already seen by the algorithm and on the distribution of the data (which does not need to be known). Still, the samples need to be identically distributed and independent or exchangeable. It is an ``online'' method using the $n$ previous samples and labels to predict the label (the training samples) and the confidence region of the next one (the test samples). For this, one calculates the non-conformity score for every new sample and decides whether the new sample lies in the prediction region or not, which is defined as

$

    \mathbb{P}(X_{n+1} \in \Gamma^{\alpha}(z_1,...,z_n,y_{n+1})) \geq  1- \alpha

$

with $\Gamma^{\alpha}(z_1,...,z_n,y_{n+1})$ denoting the prediction region and $z_{i} := (y_i,x_i)$ denoting the samples. Then, the method updates $\Gamma^{\alpha}$ at every step, also see [FZV23].


Recently, some works used conformal prediction to establish confidence intervals in a similar fashion to ours, e.g., *Conformal Prediction Masks: Visualizing Uncertainty in Medical Imaging* [K23]. They created conformal prediction-based uncertainty masks for imaging tasks. However, we see a few caveats: 1) computational costs -- The approach requires training an additional mask model and performing a calibration step, which adds computational overhead compared to methods that directly output uncertainty estimates. Instead, we don't need to perform any additional training, and the debiased method is fast. 2) Sensitivity to the choice of divergence measure -- the results might vary significantly depending on the chosen divergence measure, and it's not always clear which measure is most appropriate for a given task when one wants to quantify the uncertainty. Our debiased method has a very concrete metric on the uncertainty. 3) each value of the prediction mask is defined independently from other values. Hence, it requires the user to specify a risk level for each pixel, which is cumbersome, especially in high-dimension. Our method does not require to define a risk level for each pixel. Although it is sufficient to define a global one, our method is flexible enough to handle pixel-wise different significance levels, if desired.

Debiased methods, we believe, may be preferred when accurate coefficient estimation and hypothesis testing are primary goals, whereas conformal prediction excels in scenarios where reliable prediction intervals (e.g., images of new patients based on images of previous patients) are crucial. Overall, we believe that the two methods are no competitors, but rather, combining the advantages of the two methods (the generality of CP with the sharpness of debiased methods) is an interesting future research direction. We will add a refined version of this discussion and a broader literature review for UQ for regression problems in the final version.

[FZV23] Fontana, Matteo, Gianluca Zeni, and Simone Vantini. "Conformal prediction: a unified review of theory and new challenges." Bernoulli 29.1 (2023): 1-23.

[K23] Kutiel, Gilad, et al. "Conformal prediction masks: Visualizing uncertainty in medical imaging." International Workshop on Trustworthy Machine Learning for Healthcare. Cham: Springer Nature Switzerland, 2023.

• W2. The figures are too small, making them hard to interpret. This is compounded by the size of the text in the images. Their presentation should be rethought and fixed.
  • Thanks. We will increase the size of the figures for the final version, which is possible since we will get one page more. Please check the attached pdf with some of the new figures.

Thank you for your review and question about generalizing our result to other distributions. We carefully address it below. We would like to re-emphasize that our method is general enough and allows, for the first time, to quantify uncertainty when we do not have access to the ground truth and estimates of our solution (as in the case of complicated neural networks). We illustrate the method by using a SOTA network for inverse problems (It-Net by Genzel et al. (2022a)), but the method can be applied to many scenarios. That said, we would really appreciate it if the reviewer could raise the score to properly acknowledge the generality, rigor, and broad applicability of our method.

• Q1. Can similar to Theorem 3 results be provided for other well-known distributions? Maybe some heavy-tail dist?
  • Thanks for the excellent question. Our Theorem 3 was focused on the Gaussian case because we observed that the remainder term follows a Gaussian distribution in many MRI settings. However, Theorem 3 can indeed be generalized once the distribution of $R$ is known (even for a heavy-tailed distribution as illustrated below with a complex t-distribution), and we will include this proof and discuss how to generalize it in the final version. More precisely, we can bound the estimation error $\vert \hat{x}^u_j - x_j\vert$ analogously to the proof of Theorem 2 by $\mathbb{P}(\vert \hat{x}^u_j - x_j\vert \geq r_j(\alpha)) = \mathbb{P}( \vert W_j + R_j \vert > r_j(\alpha)) \leq \mathbb{P}( \vert W_j \vert > r_j^W(\alpha)) + \mathbb{P}( \vert R_j \vert > r_j^R(\alpha))$. The distribution of $\mathbb{P}( \vert W_j \vert > r_j^W(\alpha))$ is determined by the Gaussian noise and can be similarly computed as in Theorem 2. But if, instead, $R$ follows a heavy-tailed distribution, we can still compute $\mathbb{P}( \vert R_j \vert > r_j^R(\alpha))$ and choose $\alpha$ such that $\mathbb{P}(x_j \in C_j(\alpha)) \geq 1-\alpha$ holds for a given radius $r_j(\alpha)$. We just established the following theorem to illustrate our point for a certain heavy-tailed distribution.

Theorem Let $\hat{x}^u\in\mathbb{C}^N$ be a debiased estimator for $x\in\mathbb{C}^N$ with a remainder term following a complex t-distribution with a degree of freedom $\nu>2$, i.e. $R\sim\mathcal{C}t_\nu(0,\Sigma_R)$. Set $\eta_j:=(\Sigma_R)_{jj}$. Then, $C_j(\alpha)=\\{ z \in\mathbb{C}\mid \vert z-\hat{x}^u_j\vert \leq r_j(\alpha)\\}$ with radius

$

    r_j(\alpha) = \frac{\sigma(M\hat{\Sigma}M^\ast)_{jj}^{1/2}}{\sqrt{m}}\sqrt{\log\left(\frac{1}{\gamma_j \alpha}\right)} + \sqrt{\frac{\eta_j\nu}{2}}\sqrt{(1-\gamma_j)^{-2/\nu} \alpha^{-2\nu} - 1}.

$

is valid, i.e. $\mathbb{P}\left( x_j \in C_j(\alpha)\right)\geq 1-\alpha$.

Proof We can bound the estimation error $\vert \hat{x}^u_j - x_j\vert$ analogously to the proof of Theorem 2 by

$

    \mathbb{P}(\vert \hat{x}^u_j - x_j\vert \geq r_j(\alpha)) &= \mathbb{P}( \vert W_j + R_j \vert > r_j(\alpha)) \leq \mathbb{P}( \vert W_j \vert > r_j^W(\alpha)) + \mathbb{P}( \vert R_j \vert > r_j^R(\alpha))

$

The distribution of $\mathbb{P}( \vert W_j \vert > r_j^W(\alpha))$ is determined by the Gaussian noise and can be computed as

$

    r_j^W(\alpha) = \frac{\sigma(M\hat{\Sigma}M^*)_{jj}^{1/2}}{\sqrt{m}}\sqrt{\log\left(\frac{1}{\gamma_j\alpha}\right)}

$

similarly to Theorem 2. If $R\sim\mathcal{C}t_\nu(0,\Sigma_R)$, then the marginal distribution $R_j$ is also complex t-distributed, i.e. $R_j \sim(0,\eta_j)$, see [OTKP12]. Moreover, the probability density function of $\vert R_j\vert$ is

$

    f(r) = \frac{2r}{\eta_j}\left(1+\frac{2r^2}{\eta_j \nu}\right)^{-(\nu/2+1)}.

$

Hence,

$

\mathbb{P}( \vert R_j \vert > r_j^R(\alpha))
= \int\limits_{r_j(\alpha)}^{\infty} f(r) dr
= \left[ -\frac{1}{\left(\frac{2r^2}{\eta_j \nu}+1\right)^{\nu/2}} \right]_{r_j(\alpha)}^{\infty}
= \frac{1}{\left(\frac{2r^2}{\eta_j \nu}+1\right)^{\nu/2}}.

$

Setting this probability equal to $(1-\gamma_j)\alpha$ requires

$

r_j^R(\alpha) = \sqrt{\frac{\eta_j \nu}{2}}\sqrt{(1-\gamma_j)^{-2/\nu} \alpha^{-2\nu} - 1}.

$

Since we assumed $\nu>2$ the inequality $(1-\gamma_j)^{-2/\nu} \alpha^{-2\nu}>1$ holds. Combining $r_j(\alpha) = r_j^W(\alpha)+r_j^R(\alpha)$ concludes the proof.

[OTKP12] - Esa Ollila, David E. Tyler, Visa Koivunen, and H. Vincent Poor. Complex Elliptically Symmetric Distributions: Survey, New Results and Applications. IEEE Transactions on Signal Processing, 60(11):5597–5625, 2012

We thank you for your meticulous examination of our paper and for offering valuable feedback and criticism. We are particularly thankful for suggestions to improve clarity and for acknowledging that our work is likely to benefit a variety of high-dimensional regression applications. We address the weaknesses and questions below, and for the final version, we will make the figures larger since we also have one more page available. Please see also the attached pdf with some of the figures.

• W1. Please make the figures larger. The figures are missing axis labels and/or legends. Also, the tick and axis labels are too small. Please label individual subfigures in addition to describing them in the figure caption (e.g., "(a) w/o data adjustment" for Figure 1). For subfigures 1(d) and 1(e), and similar figures throughout the text, it would be helpful to overlay the confidence level $1-\alpha$ in a horizontal line. For Figure 3, please display (b) and (c) on the same y-axis scale.

  • Thank you for all the comments to improve the readability of our figures. The new ones are attached in the pdf. We will change all of them (and enlarge them) in the final version since we have one page more.

• W2. L49: confusing phrasing, "when the dimensions of the problem grow" to describe the asymptotic setting''.

  • With this sentence, we want to express that the dimension of the ground truth $N$, as well as the dimension of the measurements $m$, tends towards infinity with a specific ratio. This is a common assumption in the high-dimensional statistics literature. See, e.g., the papers:

    • Javanmard, A. and Montanari, A.. Hypothesis Testing in High-Dimensional Regression under the Gaussian Random Design Model: Asymptotic Theory. IEEE Trans. on Inform. Theory, 60(10):6522–6554, 2014
    • van de Geer, S., Bühlmann, P., Ritov, Y., and Dezeure, R. (2014). On asymptotically optimal confidence regions and tests for high-dimensional models. The Annals of Statistics, 42(3).

    However, in the introduction, we want to avoid explanations that are too technical. We will add a short sentence clarifying this point.

• Q1. Is the method prone to inefficiency or overcoverage of the CIs? It would be great to see some discussions in the experiments section as to where we could gain precision, for instance, from the optimization of $\gamma$, and also refer to Section A in the main text.

  • That is a great question. One main advantage of our method is its generality, i.e., that it does not require assumptions on the distribution of the ground truth data except for the existence of the second moment. We handle this general case by exploiting an empirical version of Chebyshev's inequality for the remainder term, which is sharp. Hence, our approach can deal with distributions that behave badly. If you have a nice distribution, then our method might be prone to overcoverage. This is the price we pay for its generality. One way to overcome this trade-off is to use higher moments in the estimation process, e.g., starting with the fourth moment. Regarding the choice of $\gamma$, we will discuss it and refer to Section A in the main text. In particular, we will provide numerics for different choices of $\gamma$, along with a discussion about them in the experiments section. Thank you for this great suggestion.

• Q2. What is meant by the "image support?" Could the authors please elaborate in general on what means and illustrate in the case of the MRI images for some selected?

  • Sorry, we will clarify this in the paper. The image support consists of all pixels whose value is nonzero, i.e., for $x\in\mathbb{C}^N$, the (image) support is $\{i\in [1,..., N] : x_i\neq 0\}$. Roughly speaking, the support of an image contains all non-black pixels.

• Q3. Why is in $\vert W_j\vert \sim\operatorname{Rice}$ L595 and not half-normal?

  • Indeed, for a real normal variable, the absolute value would be half-normal. However, our method allows for a more general setting (since we could measure the phase of an MR image and get body movement detection), and we assume the variable $W_j$ to be complex normal, resulting in a Rice distribution for $\vert W_j\vert$.