Matrix Denoising with Doubly Heteroscedastic Noise: Fundamental Limits and Optimal Spectral Methods

We thank Reviewer Fvn7 for the insightful comments and positive evaluation. Below we address the comment concerning our various assumptions.

We agree with the reviewer that the Gaussian-ness assumptions, the assumption that $n,d$ are of comparable order, and the assumption that the empirical spectral distributions of $\Xi, \Sigma$ converge are mild, and they basically allow the use of asymptotic methods.

We now make a few comments on the assumption that $\Xi, \Sigma$ are known, and we will add this discussion to the revision.

• In some settings, it is possible to estimate $\Xi, \Sigma$ (even consistently). This is the case if such matrices possess additional structures, e.g., if they are sparse [R4], if their inverses are sparse [R5] or if they are circulant or Toeplitz [R6]. We refer to the survey [R7] for detailed results on estimating structured high-dimensional covariance matrices.

• Recently [R8] addresses the challenge of unknown covariances by considering a modified model where one additionally observes an independent copy of noise. The statistician can then estimate the covariance from the noise-only observation and use it as a surrogate of the true covariance for estimating the signals from the spiked model. It’s possible to derive similar results in the doubly heteroskedastic setting considered in our paper.

• Finally, if the covariances are completely unknown, then our model (with Gaussian priors) is equivalent to a spiked matrix model with a certain bi-rotationally invariant noise. This problem is expected to exhibit rather different behaviors than when the covariances are known. See [R9] and [R10] for recent progress on understanding the statistical and computational limits for such models.

[R4] T. T. Cai, H. H. Zhou, "Optimal rates of convergence for sparse covariance matrix estimation", Annals of Statistics, 2012.

[R5] M. Yuan, "High dimensional inverse covariance matrix estimation via linear programming", Journal of Machine Learning Research, 2010.

[R6] T. T. Cai, Z. Ren, H. H. Zhou, "Optimal rates of convergence for estimating Toeplitz covariance matrices", Probability Theory and Related Fields, 2013.

[R7] T. T. Cai, Z. Ren, H. H. Zhou, "Estimating structured high-dimensional covariance and precision matrices: Optimal rates and adaptive estimation", Electronic Journal of Statistics, 2016.

[R8] M. Gavish, W. Leeb, E. Romanov, "Matrix denoising with partial noise statistics: optimal singular value shrinkage of spiked F-matrices", Information and Inference: A Journal of the IMA, 2023.

[R9] J. Barbier, F. Camilli, M. Mondelli, Y. Xu, "Information limits and Thouless-Anderson-Palmer equations for spiked matrix models with structured noise." arXiv preprint arXiv:2405.20993, 2024.

[R10] R. Dudeja, S. Liu, J. Ma, "Optimality of Approximate Message Passing Algorithms for Spiked Matrix Models with Rotationally Invariant Noise", arXiv preprint arXiv:2405.18081, 2024.

We thank Reviewer S4fP for the positive evaluation. Below we address the comments and questions.

Strictly speaking, you do not show that whitening fails, only that the whitened matrix does not match the proposed AMP approach (lines 262-265).

We will make the following 2 clarifications regarding lines 262-265.

1. One can repeat the analysis of an AMP that operates on the whitened matrix $\Xi^{-1/2} A \Sigma^{-1/2}$. The fixed point equations of the resulting state evolution do not match the information-theoretically optimal one in (4.5). In particular, the weak recovery threshold coming out of this approach is strictly larger than the optimal one in (4.6), as long as at least one of $\Xi, \Sigma$ is not a multiple of the identity. Since these derivations led to suboptimal results, the details were left out from the paper.
2. Our optimal spectral estimator is motivated by Bayes-AMP. If one instead considers spectral estimators associated with the whitened matrix $\Xi^{-1/2} A \Sigma^{-1/2}$, then the resulting weak recovery threshold and estimation error (overlap and matrix MSE) are all suboptimal. The exact asymptotic values of these quantities can be retrieved from Leeb–Romanov [40] and Leeb [41]. Numerical results and the corresponding theoretical predictions for whitened spectral estimators are shown in Figure 2; note its suboptimality compared to our spectral estimator and the information-theoretic limit.

To summarize, the whitening approach is provably suboptimal for both AMP and spectral estimators.

Initially you say that you believe (5.1) implies \sigma_2^{\star} < 1, and later you say that you believe these conditions are equivalent. Could you please clarify which one it is?

We apologize for the confusion, and we will clarify this point in the revision. We believe that these two conditions are actually equivalent. However, the proof of Theorem 5.1 only requires one direction: (5.1) implies $\sigma_2^* < 1$. Therefore, we formally only make the minimal conjecture of (5.1) implying $\sigma_2^* < 1$.

We thank Reviewer hDgq for carefully reading the manuscript and for the insightful comments and suggestions. Below we address each point raised in the review separately.

I believe one issue in the writing is its lack of clarity in stating which results are completely rigorous and which aren't

While our approach is indeed inspired by statistical physics, we will clarify that all our results (i.e., Proposition 4.1, Theorem 4.2, Corollary 4.3, Theorem 4.4, Theorem 5.1 and Corollary 5.2) are mathematically rigorous, with the only technical condition being “(5.1) implies $\sigma_2^* < 1$” in Theorem 5.1 that we only managed to verify numerically, but not analytically.

For the information-theoretic results, our proof uses the Gaussian interpolation method which is a rigorous technique originating in physics.

For results on spectral estimators, the analysis is inspired by Bayes-AMP. A similar approach (albeit for different learning problems) is put forward e.g. in [R1] and [R2]. However, in sharp contrast with the two works above which only provide heuristics, our results are completely rigorous. In fact, a key element of novelty in our paper is precisely to give an exact asymptotic characterization of spectral estimators via AMP tools. Thus, the heuristics on pages 8-9 are just to offer the readers an intuition on how the spectral estimators given in (5.4) arise from Bayes-AMP. We will make this clear in the revised version.

I think it would greatly improve the readability of figures 2 and 3 to have the error on the mean instead of the std.

In Figures 2-3, we report the mean averaged over 20 trials, as well as the error bars representing 1 standard deviation from the mean. Please let us know if you have any suggestions on how to improve the readability of the figures.

I would add a sentence describing more explicitly the difference and respective advantages of (3.1) (where the spectrum is O(1)) and (4.1) (where the elements of Y are O(1)).

The different scalings in (3.1) and (4.1) are purely for mathematical convenience. For the analysis of spectral estimators, it’s more convenient to have matrices whose operator norms are of constant order. For the information-theoretic analysis, it is customary in the literature to have Hamiltonians on the order of $n$ and then normalize the free energy (i.e., the expected log partition function) by $1/n$. In principle, all results and proofs can be written under a single scaling. We will clarify this in the revised version.

Small typos

Thank you for spotting the typos, we will revise accordingly.

Usage of $\overline{\Xi}, \overline{\Sigma}$

We will make it clear that all expectations involving these random variables are computed as integrals against the limiting spectral distributions of $\Xi, \Sigma$, as is common in the random matrix theory literature.

Would it be possible to look at the singular values of $A$ before and after the pre-processing? Can we clearly see a spike emerging if 5.1 is true by doing the pre-processing?

Yes, we clearly see a spike emerging if (5.1) holds by doing the pre-processing. Thank you for this excellent suggestion. We have added a plot in the PDF attached to the global response that shows the presence of spectral outliers in $A^*$, as well as the absence of such outliers in $A$. We will add this plot to the revision.

Initialising the estimator at random from the prior should allow it to have non-zero overlap in O(log(d)) steps

For some problems, AMP with random initialization may be able to attain performance similar to AMP with warm start (such as spectral initialization). Such a phenomenon has been empirically observed for various PCA and regression problems, including ours. However, proving such a behavior largely remains open. To our knowledge, the only progress is the recent work of Li, Fan and Wei [R3] for $\mathbb{Z}_2$ synchronization, i.e., rank-1 matrix estimation with Rademacher prior and GOE noise (a setting much simpler than the one considered here).

Could you run AMP after using the spectral estimator and put the additional lines and simulation dots in figures 2, 3?

All experiments in Figures 2,3 are for Gaussian priors and, in this case, our spectral estimators attain the same asymptotic performance as Bayes-AMP, as the heuristics on page 8-9 suggest. There is therefore no advantage of running AMP with spectral initialization, as opposed to the spectral method alone. Furthermore, in the special case of Figure 2a, our spectral estimators alone are information-theoretically optimal.

[R1] A. Maillard, F. Krzakala, Y. M. Lu, L. Zdeborová, "Construction of optimal spectral methods in phase retrieval", Mathematical and Scientific Machine Learning, 2022.

[R2] E. Troiani, Y. Dandi, L. Defilippis, L. Zdeborová, B. Loureiro, F. Krzakala, "Fundamental limits of weak learnability in high-dimensional multi-index models", arXiv preprint arXiv:2405.15480, 2024.

[R3] G. Li, W. Fan, Y. Wei, "Approximate message passing from random initialization with applications to Z 2 synchronization", Proceedings of the National Academy of Sciences, 2023.