1 Introduction↩︎

Covariance estimation is a classical task of critical importance across a variety of scientific and engineering fields. Applications range from machine learning, graphical models, and compressed sensing [2] to wireless communications, signal processing [3]–[7], and image analysis [8]–[10]. Furthermore, it is widely utilized in finance [11]–[16], bioinformatics, genomics, microarrays [12], [17]–[22], and generic random matrix theory [23].

The most well-known estimator is the classical sample covariance, \(\hat{\Sigma}\), which has been extensively studied over the last several decades. Typically, the primary interest lies in its deviation from the true covariance, \(\Sigma\), and how this deviation is impacted by sample complexity [1], [24]–[31]. Because these underlying problems are challenging, analytical error estimates are usually qualitative bounds focused on establishing correct scaling orders with respect to the ambient dimension, \(d\), and sample size, \(n\).

Here, we introduce a different methodology that allows for precise characterizations beyond scaling orders, enabling explicit descriptions of how sample complexity impacts estimation error. For example, in practically dominant large-dimensional contexts with a sample complexity ratio \(\alpha=\frac{n}{d}\), our results can determine the exact value of the sample covariance error’s spectral norm, \(|\hat{\Sigma}-\Sigma|_2\). This allows us to answer highly practical questions—such as how much the error decreases if the sample size is doubled or tripled—which are typically difficult to resolve within traditional scaling-order frameworks.

1.1 Relevant prior work↩︎

Classical scenarios: Various metrics of sample covariance error have been of interest. Early studies, dating back to the 1990s, focused on quantifying the error’s operator norm in the context of volume estimation and convex bodies [24]–[26]. The initial upper-bounding estimate \(O(\frac{d}{\sqrt{n}})\) [24] related to isotropic vectors was successively improved in [24]–[29] and eventually brought to the expected \(O(\sqrt{\frac{d}{n}})\) in [30] for a generic class of ensembles, including log-concave ones. Similar results of \(O(\sqrt{\frac{d}{n}}, \frac{d}{n})\) were obtained in [32] for a generic non-isotropic scenario where the true covariance matrix is different from identity (throughout this presentation, we are typically interested in sufficiently over-sampled scenarios, which are likely the most practically relevant; in such contexts, \(O(\sqrt{\frac{d}{n}}, \frac{d}{n}) = O(\sqrt{\frac{d}{n}})\)). This is scaling-wise optimal in typical scenarios where data are highly heterogeneous. However, when facing more homogeneous data where the covariances reside in a small dimensional subspace (i.e., they are of a fairly small rank), one would expect that the \(d\)-dependence might be improved.

Rank deficient \(\Sigma\): A strong effort has been made to not only improve, but completely remove the \(d\)-dependence. Relying on the so-called effective rank as a (dimensionless) measure of rank deficiency, [32], [33] observed that in the above results \(d\) can be decreased to the product of effective rank and \(\log(d)\) if the support of the underlying distribution is limited to the Euclidean ball. Analogous results were obtained through the consideration of covariance estimation with missing observations [34]. Nonetheless, while weaker, the \(d\)-dependence still remained.

The complete, formal removal of \(d\)-dependence was finally achieved in a subgaussian context in the breakthrough work [1], where a generic chaining methodology [35]–[37] (see also, [38]–[41]) was utilized to obtain an error characterization dependent solely on the effective rank and \(n\) (see also, [42] for a reproof). Subsequently, the results obtained in [1] have been significantly extended so that their correspondingly modified forms cover smooth functions of the error’s operator norm, spectral covariance projectors, and shift models [43]–[47] (for extensions in different technical directions, see [48], [49]).

The results of [50], [51] consider the so-called Frobenius norm as another type of error metric and prove effective rank bounding estimates analogous to the operator norm ones from [1]. For the importance and relevance of the Frobenius norm error metric in hypothesis testing, see, e.g., [52], [53].

Various (distributional, parametric, structural, contaminated, etc.) extensions: Many of the results discussed can be adapted to hold in similar or even improved forms when \(\Sigma\) is a priori known to possess a particular structure. These adaptations include both potential estimator restructuring and refinements in theoretical analyses. A variety of structures have been considered throughout recent literature, including: sparse [16], [54]–[56], block [57], [58], Toeplitz [59], [60], and Kronecker [61], [62]. For more on general covariance structuring and its relevance, please see [63].

Different distributional aspects, including anisotropic and heavy-tails, have also been of interest [64]–[71]. Additionally, deviations from distributional identicalness—with a particular emphasis on independent scenarios with different variance profiles—have been addressed in [72]–[75]. In this regard, particularly valuable progress has been achieved in studying matrix concentrations via free probability [76], [77] (for universality aspects, please refer to [78]).

One often encounters scenarios where the available data is gathered in a non-ideal fashion. In these situations, the robustness properties of estimators become of prevalent interest. Recent literature has addressed robust covariance estimation in many non-ideal scenarios. This includes typical cases where the data are incomplete or missing [34], [79]–[83], or corrupted and contaminated in various ways [67]–[69], [84]–[87]. Corresponding robust mean estimations are equally relevant. Many of the techniques developed for handling either covariance or mean estimation can be reutilized to address the other. For recent progress, please see, for example, [64], [88]–[93]. More information on classical corruptive and contaminated models can be found in standard robust estimation references, such as [89], [94]–[96].

There is a host of problems where data is structured differently and similar, but often more complicated, estimation tasks are faced. These include tensor or cross-covariance considerations [97]–[101], multi-inference alignment [102], [103], and (sparse) principal components/SVD analysis [104]–[106]. In each of these scenarios, one can establish analogous types of estimators and characterize the residual estimation error in a fashion similar to traditional covariance estimation.

The problems we are studying here are to some degree related to topics in statistical inference and spectral analysis as well. Characterizing the properties of spiked models, recovering low-rank matrices/tensors from noisy observations, and analyzing the spectra of deformed random matrices are prominent examples where excellent progress has been made over the last decade. Determining their key features, residual estimation errors, BBP phase transitions [77], [107]–[118], and the locations of spectral edges—along with associated large deviations principles (LDP) ([119]–[124])—typically requires a highly nontrivial transition from scaling to precise analyses.

1.2 Our contributions↩︎

As stated earlier, we pursue a different direction here, move beyond scaling estimates, and focus on the precise characterization of the error’s spectral norm. In particular, we consider a proportional large-dimensional setup with \(\alpha= \frac{n}{d}\) remaining fixed as \(d\) and \(n\) grow. For centered Gaussian vectors with true covariance \(\Sigma\), we consider the sample covariance \(\hat{\Sigma}\) and determine the limiting average spectral norm of the error \(\delta(\alpha) =\lim_{d\rightarrow\infty} \|\hat{\Sigma} - \Sigma\|_2\).

We achieve this by introducing a generic framework based on Random Duality Theory (RDT). First, a closed-form, explicit RDT-based upper bound is established (Sections 3.1–3.3). We then create a novel bilinear-quadratic RDT lower-bounding mechanism (Section 3.4). By combining this mechanism with a 2-replica systems bounding strategy, we ultimately show that it matches the upper bounds (Section 3.4.2). Numerical evaluations and simulations are also conducted, showing excellent agreement with the theoretical predictions even for relatively small problem sizes on the order of thousands (Figures 1–3 and Section 3.6).

Almost all of the results discussed in the previous section are of the scaling-order type. To move beyond these scaling barriers, a fundamentally different analytical mechanism was needed. A key benefit of this newly developed machinery is its versatility; it can be used to handle many, if not all, of the scenarios studied in prior literature and discussed in Section 1.1. Since this is an introductory paper, we chose the most classical variant of the problem. Extensions to encompass more advanced scenarios rely on the same concepts presented here, but are technically problem-specific and will be discussed elsewhere.

2 Sample covariance error – mathematical preliminaries↩︎

We start by introducing precise definitions of mathematical objects used throughout the presentation. Let \(n\in{\mathbb{N}}\) and \(d\in{\mathbb{N}}\) be two positive integers such that \[\label{eq:amat1a0} \alpha \triangleq \lim_{d\rightarrow\infty} \frac{n}{d}= const.\tag{1}\] We follow the traditional literature convention and refer to the high-dimensional regime characterized by (1 ) as linear or proportional. Let \(\bar{{\boldsymbol{x}}}\in{\mathbb{R}}^{d\times 1}\) be a \(d\)-dimensional vector comprised of centered Gaussian elements. Moreover, let the associated covariance be \(\Sigma\), i.e., let \[\label{eq:amat1a0a0} {\mathbb{E}}\bar{{\boldsymbol{x}}}\bar{{\boldsymbol{x}}}^T =\Sigma.\tag{2}\] We then also write \(\bar{{\boldsymbol{x}}}\sim {\mathcal{N}}(0,\Sigma)\). A classical way to estimate \(\Sigma\) is to collect a sample \(\{ \bar{{\boldsymbol{x}}}^{(1)},\bar{{\boldsymbol{x}}}^{(2)},\dots,\bar{{\boldsymbol{x}}}^{(n)}\}\) of \(n\) independent draws of \(\bar{{\boldsymbol{x}}}\) and to average them out. Formally, for \(\bar{{\boldsymbol{x}}}^{(i)}\sim {\mathcal{N}}(0,\Sigma),i\in\{1,2,\dots,n\}\), one sets \[\label{eq:amat1a0a1} X \triangleq \begin{bmatrix} \left (\bar{{\boldsymbol{x}}}^{(1)} \right )^T \\ \left (\bar{{\boldsymbol{x}}}^{(2)} \right )^T \\ \vdots \\ \left (\bar{{\boldsymbol{x}}}^{(n)} \right )^T \end{bmatrix},\tag{3}\] and considers the sample covariance matrix \(\hat{\Sigma} \triangleq \frac{1}{n}X^TX\) as an approximation to \(\Sigma\). In other words, one hopes that \[\label{eq:amat1a0a2} \hat{\Sigma} \triangleq \frac{1}{n} X^TX \approx \Sigma \quador\quad \|\hat{\Sigma} - \Sigma\| = \left \| \frac{1}{n}X^TX - \Sigma \right \|is small ,\tag{4}\] where \(\|\cdot\|\) is the norm of choice. Naturally, as the sample size, \(n\), grows (which also means that \(\alpha\) increases), the quality of the approximation is expected to improve. Precisely characterizing this improvement remains one of the fundamental challenges in statistical estimation and probability theory.

A major breakthrough in this area was achieved in [1], where the error’s operator norm was shown to behave in the following way \[\label{eq:amat1a0a3} \delta_n \triangleq {\mathbb{E}}\|\hat{\Sigma} - \Sigma\| \sim \|\Sigma\| \max \left (\sqrt{\frac{r(\Sigma)}{n}} , \frac{r(\Sigma)}{n} \right ),\tag{5}\] where \[\label{eq:amat1a0a4} r(\Sigma) \triangleq \frac{\left ({\mathbb{E}}\|\bar{{\boldsymbol{x}}}\|\right )^2}{\|\Sigma\|}.\tag{6}\] Accompanying concentration inequalities are established in [1] as well. The estimation error introduced in (5 ) is a function of \(n\), and larger sample sizes expectedly produce smaller errors (i.e., as \(n\) increases, \(\delta\) decreases). The characterizations in (5 ) and (6 ) provide the correct error scaling, highlighting the roles of the sample size, \(n\), and the norm of the true covariance, \(\|\Sigma\|\). While this provides foundational information regarding the quality of \(\hat{\Sigma}\), practical scenarios often demand more. Specifically, when gathering abundant data is costly or difficult, precisely understanding the tradeoff between sample size and accuracy becomes a necessity.

In the remainder of the paper, we address this demand and determine the underlying tradeoff. In particular, specializing for concreteness to the typical spectral norm choice, we precisely characterize \(\delta(\alpha)\), defined as \[\label{eq:amat1a0a5} \delta(\alpha)\triangleq \lim_{n\rightarrow \infty }\delta_n = \lim_{n\rightarrow \infty } {\mathbb{E}}\|\hat{\Sigma} - \Sigma\|_2 = \lim_{d\rightarrow \infty } {\mathbb{E}}\|\hat{\Sigma} - \Sigma\|_2,\tag{7}\] where the last equality holds based on (1 ).

Let the output of function \(\lambda_i(\cdot)\) be the \(i\)-th smallest eigenvalue of its argument. Then one has \[\label{eq:amat1a0a6} \|\hat{\Sigma} - \Sigma\|_2 = \max(\lambda_n(\hat{\Sigma} - \Sigma),|\lambda_1(\hat{\Sigma} - \Sigma)|).\tag{8}\] For the time being, we focus on the largest eigenvalue of \(\left (\hat{\Sigma} -\Sigma \right )\), i.e., on \(\lambda_n\left (\hat{\Sigma} -\Sigma \right )\) (later on, in Section 3.5, we will see that such an approach suffices to determine \(\delta(\alpha)\)). To that end, we denote the unit sphere in \({\mathbb{R}}^d\) by \({\mathbb{S}}^d=\{{\boldsymbol{x}}|\|{\boldsymbol{x}}\|_2=1\}\) and write \[\begin{align} \label{eq:inteq1ad0} \lambda_n\left (\hat{\Sigma} -\Sigma \right ) = \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d}{\boldsymbol{x}}^T \left (\hat{\Sigma} -\Sigma \right ){\boldsymbol{x}} = \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d}{\boldsymbol{x}}^T \left (\frac{1}{n}X^TX -\Sigma \right ){\boldsymbol{x}} = \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d} \left (\frac{1}{n}{\boldsymbol{x}}^T X^TX {\boldsymbol{x}}- {\boldsymbol{x}}^T\Sigma {\boldsymbol{x}}\right ). \end{align}\tag{9}\] By definition, \(\Sigma\) is positive semi-definite (\(\Sigma\succeq 0\)) and admits the following eigen-decomposition \[\begin{align} \label{eq:inteq1ad1} \Sigma = USS^TU^T. \end{align}\tag{10}\] In (10 ), \(S\in{\mathbb{R}}^{d\times d}\) is a \(d\times d\) diagonal matrix (with roots of \(\Sigma\)’s eigenvalues on the main diagonal) and \(U\in{\mathbb{R}}^{d\times d}\) is a \(d\times d\) orthogonal matrix with \(UU^T=U^TU=I\) (to avoid sidetracking the flow of the presentation with constant mentions of special cases, we assume that the eigenvalues of \(S\) are positive and belong to an interval independent of \(d\).). Let \(A\in{\mathbb{R}}^{n\times d}\) be an \(n\times d\) matrix comprised of independent standard normal elements. We can then write the following statistical equivalent to (9 ) \[\begin{align} \label{eq:inteq1ad2} \lambda_n\left (\hat{\Sigma} -\Sigma \right ) & = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d} \left (\frac{1}{n}{\boldsymbol{x}}^T X^TX {\boldsymbol{x}}- {\boldsymbol{x}}^T\Sigma {\boldsymbol{x}}\right ) \nonumber \\ & = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d} \left (\frac{1}{n}{\boldsymbol{x}}^T US^TU^TA^TAUSU^T {\boldsymbol{x}}- {\boldsymbol{x}}^T USS^TU^T {\boldsymbol{x}}\right ) \nonumber \\ & = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d} \left (\frac{1}{n}{\boldsymbol{x}}^T S^TU^TA^TAUS {\boldsymbol{x}}- {\boldsymbol{x}}^T SS^T {\boldsymbol{x}}\right ) \nonumber \\ & = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d} \left (\frac{1}{n}{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}- {\boldsymbol{x}}^T SS^T {\boldsymbol{x}}\right ). \end{align}\tag{11}\] Since \({\boldsymbol{x}}\in{\mathbb{S}}^d\), the third equality follows through the cosmetic change of variables, \(U^T{\boldsymbol{x}}\rightarrow {\boldsymbol{x}}\). Similarly, the fourth equality is enabled by the rotational invariance of \(A\), i.e., by the fact that \(A\) and \(AU\) have the same distribution. For a real scalar \(c\), we find it useful to set \[\begin{align} \label{eq:inteq1a} \xi(c) & = & \frac{1}{\sqrt{n}} \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c } \sqrt{{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}}. \end{align}\tag{12}\] This allows to rewrite (11 ) as \[\begin{align} \label{eq:inteq1ab00} \lambda_n\left (\hat{\Sigma} -\Sigma \right ) & = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d} \left (\frac{1}{n}{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}- {\boldsymbol{x}}^T SS^T {\boldsymbol{x}}\right ) \nonumber \\ & = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c,c} \left (\frac{1}{n}{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}- {\boldsymbol{x}}^T SS^T {\boldsymbol{x}}\right ) \nonumber \\ & = & \max_{c} \left (\max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c} \frac{1}{n}{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}- c^2\right ) \nonumber \\ & = & \max_{c} \left (\left (\xi(c)\right )^2 - c^2 \right ). \end{align}\tag{13}\] To characterize \(\xi(c)\) and ultimately \(\lambda_n\left (\hat{\Sigma} -\Sigma \right )\), we utilize Random duality theory (RDT). This is done in the following section.

3 Characterizing \(\xi(c)\) via RDT↩︎

We first recall the four key RDT principles [125], [126] (for more on further upgrades and associated algorithmic implications, see, e.g., [127], [128]):

1. Finding underlying optimization algebraic representation

2. Determining random dual

3. Handling random dual

4. Double-checking strong random duality.

We below discuss each of these principles and how they relate to the problem of our interest here.

3.1 Finding underlying optimization algebraic representation↩︎

We start by observing that \(\xi(c)\) can be rewritten in the following RDT favorable fashion \[\begin{align} \label{eq:inteq1ab0} \xi(c) & = & \frac{1}{\sqrt{n}} \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c } \sqrt{{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}} \nonumber \\ & = & \frac{1}{\sqrt{n}} \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d ,\|S{\boldsymbol{x}}\|_2=c} \| AS{\boldsymbol{x}}\|_2 \nonumber \\ & = & \frac{1}{\sqrt{n}} \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c,{\boldsymbol{y}}\in{\mathbb{S}}^n} {\boldsymbol{y}}^TAS{\boldsymbol{x}}. \end{align}\tag{14}\] Within the RDT the above represents the so-called random primal.

3.2 Determining random dual↩︎

The following theorem establishes the so-called random dual.

We also note that the first two main RDT principles were rewritten in a generic processes comparisons context in [49] and the sample covariance error with the bound from (20 ) was considered as a particular application.

3.3 Handling random dual↩︎

To handle the above random dual, we closely follow the paths traced in [125], [133], [134]. Before proceeding with the details, we find it useful to note the role of the random dual with the overall RDT mosaic. Namely, a combination of (13 ), (?? ), and (?? ) together with concentrations gives \[\begin{align} \label{eq:hrd0} {\mathbb{E}}\lambda_n\left (\hat{\Sigma} -\Sigma \right ) & = & \max_{c} \left (\left ({\mathbb{E}}\xi \right )^2 - c^2 \right ) \nonumber \\ & \leq & \max_{c} \left (\left (c + {\mathbb{E}}\frac{1}{\sqrt{n}} \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c} \left ({\boldsymbol{g}}^{(2)}\right )^TS{\boldsymbol{x}}\right )^2 - c^2 \right ) \nonumber \\ & = & \max_{c} \left (\left (c + \frac{1}{\sqrt{n}} {\mathbb{E}}L(c) \right )^2 - c^2 \right ). \end{align}\tag{21}\] Clearly, \(L(c)\) is the key object of interest and we below study it in detail. First, we note \[\begin{align} \label{eq:hrd1} L(c) & \triangleq & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c} \left ({\boldsymbol{g}}^{(2)}\right )^TS{\boldsymbol{x}} = \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c} \sum_{i=1}^{d}{\boldsymbol{g}}_i^{(2)} {\boldsymbol{s}}_i{\boldsymbol{x}}_i = -\min_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c} \sum_{i=1}^{d}{\boldsymbol{g}}_i^{(2)} {\boldsymbol{s}}_i{\boldsymbol{x}}_i . \end{align}\tag{22}\] Then we have for the Lagrangian \[\begin{align} \label{eq:hrd2} {\mathcal{L}}= \sum_{i=1}^{d}{\boldsymbol{g}}_i^{(2)} {\boldsymbol{s}}_i{\boldsymbol{x}}_i + \gamma \sum_{i=1}^{d}{\boldsymbol{x}}_i^2 -\gamma + \gamma_0 \sum_{i=1}^{d}{\boldsymbol{s}}_i^2{\boldsymbol{x}}_i^2 -\gamma_0 c^2. \end{align}\tag{23}\] Combining (22 ) and (23 ) and utilizing strong duality we further write \[\begin{align} \label{eq:hrd3} L(c) = -\min_{{\boldsymbol{x}}}\max_{\gamma,\gamma_0} {\mathcal{L}}= -\max_{\gamma,\gamma_0} \min_{{\boldsymbol{x}}} {\mathcal{L}}. \end{align}\tag{24}\] To optimize over \({\boldsymbol{x}}\) we first find derivatives \[\begin{align} \label{eq:hrd4} \frac{d{\mathcal{L}}}{d{\boldsymbol{x}}_i} = {\boldsymbol{g}}_i^{(2)} {\boldsymbol{s}}_i + 2\gamma {\boldsymbol{x}}_i + 2\gamma_0 {\boldsymbol{s}}_i^2{\boldsymbol{x}}_i. \end{align}\tag{25}\] Equalling the above derivatives to zero gives \[\begin{align} \label{eq:hrd5} {\boldsymbol{x}}_i = -\frac{{\boldsymbol{g}}_i^{(2)} {\boldsymbol{s}}_i}{2(\gamma + \gamma_0 {\boldsymbol{s}}_i^2)}. \end{align}\tag{26}\] Plugging this back in (23 ), one finds \[\begin{align} \label{eq:hrd6} \min_{{\boldsymbol{x}}} {\mathcal{L}}= -\frac{1}{4}\sum_{i=1}^{d}\frac{\left ({\boldsymbol{g}}_i^{(2)}\right )^2 {\boldsymbol{s}}_i^2}{\gamma + \gamma_0 {\boldsymbol{s}}_i^2} -\gamma -\gamma_0 c^2. \end{align}\tag{27}\] A combination of (24 ) and (27 ) gives \[\begin{align} \label{eq:hrd7} L(c) & = & -\max_{\gamma,\gamma_0} \min_{{\boldsymbol{x}}} {\mathcal{L}} \nonumber \\ & = & -\max_{\gamma,\gamma_0} \left (-\frac{1}{4}\sum_{i=1}^{d}\frac{\left ({\boldsymbol{g}}_i^{(2)}\right )^2 {\boldsymbol{s}}_i^2}{\gamma + \gamma_0 {\boldsymbol{s}}_i^2} -\gamma -\gamma_0 c^2 \right ) \nonumber \\ & = & \min_{\gamma,\gamma_0} \left (\frac{1}{4}\sum_{i=1}^{d}\frac{\left ({\boldsymbol{g}}_i^{(2)}\right )^2 {\boldsymbol{s}}_i^2}{\gamma + \gamma_0 {\boldsymbol{s}}_i^2} +\gamma +\gamma_0 c^2 \right ) \nonumber \\ & = & \min_{\gamma_x,\gamma_0} \left (\frac{1}{4\gamma_0}\sum_{i=1}^{d}\frac{\left ({\boldsymbol{g}}_i^{(2)}\right )^2 {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} + \gamma_0(\gamma_x + c^2) \right ) \nonumber \\ & = & \min_{\gamma_x} \sqrt{\sum_{i=1}^{d}\frac{\left ({\boldsymbol{g}}_i^{(2)}\right )^2 {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} (\gamma_x + c^2) } \end{align}\tag{28}\] Law of large numbers and concentrations gives \[\begin{align} \label{eq:hrd8} \lim_{d\rightarrow\infty} \frac{1}{\sqrt{d}} {\mathbb{E}}L(c) & = & \lim_{d\rightarrow\infty} \sqrt{\min_{\gamma_x} \bar{L}}, \end{align}\tag{29}\] where \[\begin{align} \label{eq:hrd9} \bar{L} = \left (\frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} (\gamma_x + c^2) \right ). \end{align}\tag{30}\] Combining (21 ), (29 ), and (30 ), we find \[\begin{align} \label{eq:hrd10} \lim_{d\rightarrow \infty}{\mathbb{E}}\lambda_n\left (\hat{\Sigma} -\Sigma \right ) & \leq & \max_{c} \left (\left (c + \lim_{d\rightarrow \infty } \frac{1}{\sqrt{n}}{\mathbb{E}}L(c) \right )^2 - c^2 \right ) = \lim_{d\rightarrow\infty} \max_{c}\min_{\gamma_x} \delta_u(\alpha) , \end{align}\tag{31}\] where \[\begin{align} \label{eq:hrd11} \delta_u(\alpha) \triangleq \left (c + \frac{1}{\sqrt{\alpha}}\sqrt{ \bar{L}} \right )^2 - c^2 = 2 \frac{c\sqrt{\bar{L}}}{\sqrt{\alpha}} + \frac{\bar{L}}{\alpha} . \end{align}\tag{32}\] Computing \(\gamma_x\) derivative gives \[\begin{align} \label{eq:hrd12} \frac{d\bar{L}}{d\gamma_x} = -\frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{(\gamma_x + {\boldsymbol{s}}_i^2)^2} (\gamma_x + c^2) + \frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} = \frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2( {\boldsymbol{s}}_i^2 -c^2)}{(\gamma_x + {\boldsymbol{s}}_i^2)^2} . \end{align}\tag{33}\] Computing \(c\) derivative gives \[\begin{align} \label{eq:hrd13} \frac{d\delta_u(\alpha)}{d c} & = & \frac{2\sqrt{\bar{L}}}{\sqrt{\alpha}} + \frac{c}{\sqrt{\alpha \bar{L}}} \left (\frac{2c}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} \right ) + \frac{2}{\alpha d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} \nonumber \\ & = & \frac{2}{\sqrt{\alpha}} \left ( \frac{\bar{L}}{\sqrt{\bar{L}}} + \frac{c}{\sqrt{ \bar{L}}} \left (\frac{c}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} \right ) + \frac{c}{\sqrt{\alpha} d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} \right ) \nonumber \\ & = & \frac{2}{\sqrt{\alpha}} \left ( \frac{1}{\sqrt{\bar{L}}} \left (\frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} (\gamma_x + c^2) \right )+ \frac{c}{\sqrt{ \bar{L}}} \left (\frac{c}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} \right ) + \frac{c}{\sqrt{\alpha} d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} \right ) \nonumber \\ & = & \frac{2}{\sqrt{\alpha}}\frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} \left ( \frac{\gamma_x + c^2 }{\sqrt{\bar{L}}} + \frac{c^2}{\sqrt{ \bar{L}}} + \frac{c}{\sqrt{\alpha}} \right ) \nonumber \\ & = & \frac{2}{\sqrt{\alpha}}\frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} \left ( \frac{\gamma_x + 2c^2 }{\sqrt{\bar{L}}} + \frac{c}{\sqrt{\alpha}} \right ). \end{align}\tag{34}\] Equalling the above derivative to zero gives \[\begin{align} \label{eq:hrd14} \left ( \frac{\gamma_x + 2c^2 }{\sqrt{\bar{L}}} + \frac{c}{\sqrt{\alpha}} \right )= 0 \quad \quad \Longrightarrow \quad \quad \alpha\left (\gamma_x + 2c^2 \right )^2 = c^2\bar{L} = c^2 z (\gamma_x + c^2) , \end{align}\tag{35}\] where \[\begin{align} \label{eq:hrd15} z = \frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} . \end{align}\tag{36}\] Keeping in mind \(c\geq 0\), we also observe that equalling the last term in (34 ) to zero, additionally implies \[\begin{align} \label{eq:hrd15a0} \left ( \frac{\gamma_x + 2c^2 }{\sqrt{\bar{L}}} + \frac{c}{\sqrt{\alpha}} \right )= 0 \quad \quad \Longrightarrow \quad \quad \gamma_x\leq 0 \quad and\gamma_x+ 2 c^2\leq 0 . \end{align}\tag{37}\] One then finds the following equivalent to (35 ) \[\begin{align} \label{eq:hrd16} c^4 +c^2\gamma_x -\frac{\gamma_x^2\alpha}{z-4\alpha} = 0. \end{align}\tag{38}\] Solving for \(c^2\) gives \[\begin{align} \label{eq:hrd17} c^2 = \frac{-\gamma_x \pm \sqrt{ \gamma_x^2 +4\frac{\gamma_x^2\alpha}{z-4\alpha} }}{2} = \gamma_x\frac{- 1 \pm sign(\gamma_x) \sqrt{ 1 +4\frac{\alpha}{z-4\alpha} }}{2} = \gamma_x\frac{- 1 \pm sign(\gamma_x) \sqrt{ \frac{z}{z-4\alpha} }}{2}. \end{align}\tag{39}\] Taking into account (37 ), we then have particular choice of signs that gives \[\begin{align} \label{eq:hrd17a0} c^2 = \gamma_x\frac{- 1 + \sqrt{ \frac{z}{z-4\alpha} }}{2}. \end{align}\tag{40}\]

We first set \[\begin{align} \label{eq:hrd17a1} \phi_1(\gamma_x) \triangleq \frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^4}{\left (\gamma_x + {\boldsymbol{s}}_i^2 \right )^2}, \quad \quad \phi_2(\gamma_x) \triangleq \frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\left (\gamma_x + {\boldsymbol{s}}_i^2\right )^2}, \end{align}\tag{41}\] and find \[\begin{align} \label{eq:hrd18} z & = & \frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} = \gamma_x\phi_2(\gamma_x) + \phi_1(\gamma_x) \nonumber \\ \bar{L} & = & \frac{1}{d}\sum_{i=1}^{d}\frac{ {\boldsymbol{s}}_i^2}{\gamma_x + {\boldsymbol{s}}_i^2} (\gamma_x + c^2) = z(\gamma_x + c^2) = (\gamma_x\phi_2(\gamma_x) + \phi_1(\gamma_x))(\gamma_x + c^2). \end{align}\tag{42}\] After additionally setting \[\begin{align} \label{eq:hrd19} \phi_3(\gamma_x) = \gamma_x\frac{- 1 + \sqrt{ \frac{z}{z-4\alpha} }}{2} = \gamma_x\frac{- 1 + \sqrt{ \frac{\gamma_x\phi_2(\gamma_x) + \phi_1(\gamma_x)}{\gamma_x\phi_2(\gamma_x) + \phi_1(\gamma_x)-4\alpha} }}{2}, \end{align}\tag{43}\] the following theorem summarize handling random dual.

Remark 3 suggests that the above analysis and the resulting bounds given in Theorem 3 might be tight. We discuss this in more detail next.

3.4 Double-checking strong random duality↩︎

Theorem 1 established the random dual which upper-bounds \({\mathbb{E}}\xi(c)\). Below we discuss complementary lower bounds.

3.4.1 Bilinear-quadratic mechanism↩︎

The following theorem establishes complementary random dual which lower-bounds \({\mathbb{E}}\xi(c)\). It relies on a bilinear-quadratic comparative mechanism and is the key component that enables the whole machinery developed in the paper to work.

3.4.2 Handling \({\mathbb{E}}B(c)\)↩︎

We split handling \({\mathbb{E}}B(c)\) into two parts. \({\mathbb{E}}B(c)\)’s upper-bound is discussed in the first part and its matching lower bound in the second part.

To upper-bound \(B(c)\) we closely follow procedures developed in previous sections. The following theorem summarizes main results.

As stated earlier, Theorem 2 is a special case of concepts discussed in Corollary 3 in [131] and in Corollary 4 in [132]). The machinery developed there ensures that the upper-bounding mechanism of the previous subsection is tight provided that two nontrivially overlapped (\(q\neq1\)) replicated systems cannot double the maximal value of a single system. The very same principle was utilized in the single-partite system in [135], [136] which we consider here (albeit in a more general spin configuration).

Below, we check whether this condition indeed holds. Following [131], [132], [135], for a real scalar \(t\in [0,1]\), we first consider interpolated system \[\begin{align} \label{eq:lbstr1}\boldsymbol{1-rep:}\quad\quad D(c;t)= \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c} \left (\sqrt{t}{\boldsymbol{x}}^TS^TG^{(2)}S{\boldsymbol{x}}+\sqrt{1-t}\sqrt{2}c \left ({\boldsymbol{g}}^{(2)}\right )^TS{\boldsymbol{x}}\right ). \end{align}\tag{62}\] Clearly, \(D(c;t)\) continuously interpolates between \(B(c)\) and \(\sqrt{2}L(c)\). In particular, one has \[\begin{align} \label{eq:lbstr2} D(c,1) = B(c) \quad\quadand\quad\quad D(c;0)= \sqrt{2}cL(c). \end{align}\tag{63}\] Moreover, Theorem 5 gives \[\begin{align} \label{eq:lbstr2a0} \frac{1}{\sqrt{2d}c}{\mathbb{E}}B(c) = \frac{1}{\sqrt{2d}c} {\mathbb{E}}D(c,1) \leq \frac{1}{\sqrt{2d}c} {\mathbb{E}}D(c,t) \leq \frac{1}{\sqrt{2d}c} {\mathbb{E}}D(c;0)= \frac{1}{\sqrt{d}} {\mathbb{E}}L(c). \end{align}\tag{64}\] Then consider a set of replica pairs \[\begin{align} \label{eq:lbstr3} \bar{{\mathcal{X}}} = \left \{ ({\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)})|{\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)}\in{\mathbb{S}}^d, \|S{\boldsymbol{x}}^{(1)} \|_2= \|S{\boldsymbol{x}}^{(2)} \|_2=c,\left ({\boldsymbol{x}}^{(1)}\right )^TS^TS{\boldsymbol{x}}^{(2)}=qc^2 \right \} , \end{align}\tag{65}\] and associate to it the following \(q\)-overlapped replicated system \[\begin{align} \label{nkcawysi}\boldsymbol{2-q-rep:}\quad\quad D^{(2)}(c;t)= \max_{({\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)})\in\bar{{\mathcal{X}}}} \sum_{i=1}^{2} \left (\left ({\boldsymbol{x}}^{(i)}\right )^TS^TG^{(2)}S{\boldsymbol{x}}^{(i)} +\sqrt{1-t}\sqrt{2}c \left ({\boldsymbol{g}}^{(2)}\right )^TS{\boldsymbol{x}}^{(i)} \right ). \end{align}\tag{66}\] The above principle – nontrivially overlapped replicated system cannot double the free energy – means that for any \(t\in[0,1)\) and \(q\in(-1,1)\) \[\begin{align} \label{eq:lbstr4}\boldsymbol{2-q-rep}<\mathbf{2}\times(\boldsymbol{1-rep}). \end{align}\tag{67}\] In mathematical terminology one then has for any \(t\in[0,1)\) \[\begin{align} \label{eq:lbstr5}\boldsymbol{2-q-rep}<\mathbf{2}\times(\boldsymbol{1-rep})\quad\quad \Longleftrightarrow \quad\quad \min_{q\in(-1,1)} \left ( \lim_{d\rightarrow \infty} \frac{2}{\sqrt{d}} {\mathbb{E}}D(c;0) - \lim_{d\rightarrow \infty} \frac{1}{\sqrt{d}} {\mathbb{E}}D^{(2)}(c;t) \right )>0. \end{align}\tag{68}\] If this principle is indeed in place, i.e., if one can show that (68 ) holds then complete analogues to Theorem 2.4 in [136] and Theorem 5.2 in [135] are established and the remaining portions of the machineries of [135], [136] ensure \({\mathbb{E}}D(c,1)={\mathbb{E}}D(c,0)\). In other words, \[\label{eq:lbstr6} \min_{q\in(-1,1)} \left ( \lim_{d\rightarrow \infty} \frac{2}{\sqrt{2d}c} {\mathbb{E}}D(c;0) - \lim_{d\rightarrow \infty} \frac{1}{\sqrt{2d}c} {\mathbb{E}}D^{(2)}(c;t) \right )>0 \Longleftrightarrow\lim_{d\rightarrow \infty} \frac{1}{\sqrt{2d}c} {\mathbb{E}}D(c,1)= \lim_{d\rightarrow \infty} \frac{1}{\sqrt{2d}c} {\mathbb{E}}D(c,0).\tag{69}\]

To establish $ < () \(, we first prove the following theorem that upper-bounds\)D^{(2)}(c;t)$.

Keeping in mind (63 ), (64 ), and (?? ), condition in (68 ) and (69 ) will be met if for any \(t<1\) \[\begin{align} \label{eq:lbstrmt5a1a1} \lim_{d\rightarrow \infty} \frac{1}{\sqrt{d}} {\mathbb{E}}L^{(2)}(c;t) < 2 \lim_{d\rightarrow \infty}\frac{1}{\sqrt{d}} {\mathbb{E}}L(c). \end{align}\tag{76}\] We provide a characterization of \(\lim_{d\rightarrow \infty} \frac{1}{\sqrt{d}} {\mathbb{E}}L^{(2)}(c;t)\) in the following subsection.

Clearly, \(L(c)\) is the key object of interest. Below we study it in detail. First, from (?? ) and (?? ) \[\label{eq:l2lbstrthm1eq1} L^{(2)}(c;t)= \max_{({\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)})\in\bar{{\mathcal{X}}}} \left (\sqrt{t} \left ({\boldsymbol{g}}^{(1,1)}\right )^TS{\boldsymbol{x}}^{(1)} +\sqrt{t} \left ({\boldsymbol{g}}^{(1,3)}\right )^TS{\boldsymbol{x}}^{(2)} +\sqrt{1-t} \left ({\boldsymbol{g}}^{(2)}\right )^TS \left ({\boldsymbol{x}}^{(1)}+ {\boldsymbol{x}}^{(2)} \right )\right ).\tag{77}\] Let \({\boldsymbol{g}}^{(x,1)}\in{\mathbb{R}}^{d\times 1}\) have independent standard normal components. Also, let \({\boldsymbol{g}}^{(x,2)}\in{\mathbb{R}}^{d\times 1}\) have independent standard normal components. Additionally, let \[\label{eq:l2lbstrthm1eq1a0} {\mathbb{E}}{\boldsymbol{g}}_i^{(x,1)}{\boldsymbol{g}}_i^{(x,2)} = 1-t + qt \triangleq a, 1\leq i\leq d.\tag{78}\] One can then replace (?? ) with its statistical equivalent \[\label{eq:l2lbstrthm1eq1a1} L^{(2)}(c;t)= \max_{({\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)})\in\bar{{\mathcal{X}}}} \left (\left ({\boldsymbol{g}}^{(x,1)}\right )^TS{\boldsymbol{x}}^{(1)} + \left ({\boldsymbol{g}}^{(x,2)}\right )^TS{\boldsymbol{x}}^{(2)} \right ).\tag{79}\] Basic algebraic transformations allow to rewrite (79 ) as \[\begin{align} \label{eq:l2hrd1} L^{(2)}(c;t)= - \min_{({\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)})\in\bar{{\mathcal{X}}}} \left (\left ({\boldsymbol{g}}^{(x,1)}\right )^TS{\boldsymbol{x}}^{(1)} + \left ({\boldsymbol{g}}^{(x,2)}\right )^TS{\boldsymbol{x}}^{(2)} \right ). \end{align}\tag{80}\] We then have for the Lagrangian \[\begin{align} \label{eq:l2hrd2} {\mathcal{L}}& = & \sum_{i=1}^{d} \left ({\boldsymbol{g}}_i^{(x,1)}{\boldsymbol{s}}_i{\boldsymbol{x}}_i^{(1)} + {\boldsymbol{g}}_i^{(x,2)}{\boldsymbol{s}}_i{\boldsymbol{x}}_i^{(2)} \right ) + \gamma_1 \sum_{i=1}^{d}\left ({\boldsymbol{x}}_i^{(1)}\right )^2 -\gamma_1 + \gamma_2 \sum_{i=1}^{d}\left ({\boldsymbol{x}}_i^{(2)}\right )^2 -\gamma_2 \nonumber \\ & & + \gamma_{01} \sum_{i=1}^{d}{\boldsymbol{s}}_i^2\left ({\boldsymbol{x}}_i^{(2)}\right )^2 -\gamma_{01} c^2 + \gamma_{02} \sum_{i=1}^{d}{\boldsymbol{s}}_i^2 \left ({\boldsymbol{x}}_i^{(2)}\right )^2 -\gamma_{02} c^2 + \nu \sum_{i=1}^{d}{\boldsymbol{x}}_i^{(1)}{\boldsymbol{s}}_i^2{\boldsymbol{x}}_i^{(2)} -\nu q c^2. \end{align}\tag{81}\] A combination of (80 ) and (81 ) together with duality gives \[\begin{align} \label{eq:l2hrd3} L^{(2)}(c;t) = -\min_{{\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)}}\max_{\gamma_1,\gamma_2,\gamma_{01},\gamma_{02},\nu} {\mathcal{L}} \leq -\max_{\gamma_1,\gamma_2,\gamma_{01},\gamma_{02},\nu} \min_{{\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)}} {\mathcal{L}}. \end{align}\tag{82}\] To optimize over \({\boldsymbol{x}}^{(1)}\) and \({\boldsymbol{x}}^{(2)}\), we first find derivatives \[\begin{align} \label{eq:l2hrd4} \frac{d{\mathcal{L}}}{d{\boldsymbol{x}}_i^{(1)}} & = & {\boldsymbol{g}}_i^{(x,1)} {\boldsymbol{s}}_i + 2\gamma_1 {\boldsymbol{x}}_i^{(1)} + 2\gamma_{01} {\boldsymbol{s}}_i^2{\boldsymbol{x}}_i^{(1)} +\nu{\boldsymbol{s}}_i^2{\boldsymbol{x}}_i^{(2)} \nonumber \\ \frac{d{\mathcal{L}}}{d{\boldsymbol{x}}_i^{(2)}} & = & {\boldsymbol{g}}_i^{(x,2)} {\boldsymbol{s}}_i + 2\gamma_2 {\boldsymbol{x}}_i^{(2)} + 2\gamma_{02} {\boldsymbol{s}}_i^2{\boldsymbol{x}}_i^{(2)} +\nu{\boldsymbol{s}}_i^2{\boldsymbol{x}}_i^{(1)}. \end{align}\tag{83}\] As \(\gamma\)’s are scalars and problem is completely symmetric in \({\boldsymbol{x}}^{(1)}\) and \({\boldsymbol{x}}^{(2)}\), we have \(\gamma_1=\gamma_2=\gamma\) and \(\gamma_{01}=\gamma_{02}=\gamma_0\). Also, we set \[\begin{align} \label{eq:l2hrd4a0} {\boldsymbol{z}}_{1,i} & = & {\boldsymbol{g}}_i^{(x,1)} {\boldsymbol{s}}_i \nonumber \\ {\boldsymbol{z}}_{2.i} & = & {\boldsymbol{g}}_i^{(x,2)} {\boldsymbol{s}}_i \nonumber \\ {\boldsymbol{v}}_i & = & 2(\gamma + \gamma_0 {\boldsymbol{s}}_i^2) \nonumber \\ \nu_i & = & \nu {\boldsymbol{s}}_i^2. \end{align}\tag{84}\] Keeping (84 ) in mind and equalling derivatives in (83 ) to zero gives \[\begin{align} \label{eq:l2hrd4a1} {\boldsymbol{z}}_{1,i} + \nu_i {\boldsymbol{x}}_i^{(2)} + {\boldsymbol{x}}_i^{(1)} {\boldsymbol{v}}_i & = & 0 \nonumber \\ {\boldsymbol{z}}_{2,i} + \nu_i {\boldsymbol{x}}_i^{(1)} + {\boldsymbol{x}}_i^{(2)} {\boldsymbol{v}}_i & = & 0 . \end{align}\tag{85}\] Summing and subtracting the above two equations gives \[\begin{align} \label{eq:l2hrd4a2} {\boldsymbol{x}}_i^{(1)}+{\boldsymbol{x}}_i^{(2)} & = & -\frac{{\boldsymbol{z}}_{1,i}+{\boldsymbol{z}}_{2,i}}{\nu_i+{\boldsymbol{v}}_i} \nonumber \\ {\boldsymbol{x}}_i^{(1)}-{\boldsymbol{x}}_i^{(2)} & = & -\frac{-{\boldsymbol{z}}_{1,i}+{\boldsymbol{z}}_{2,i}}{\nu_i-{\boldsymbol{v}}_i}, \end{align}\tag{86}\] and \[\begin{align} \label{eq:l2hrd4a3} {\boldsymbol{x}}_i^{(1)} & = & \frac{ {\boldsymbol{z}}_{1,i} {\boldsymbol{v}}_i - {\boldsymbol{z}}_{2,i}\nu_i}{\nu_i^2-{\boldsymbol{v}}_i^2} \nonumber \\ {\boldsymbol{x}}_i^{(2)} & = & \frac{ -{\boldsymbol{z}}_{1,i} {\boldsymbol{v}}_i + {\boldsymbol{z}}_{2,i}\nu_i}{\nu_i^2-{\boldsymbol{v}}_i^2} . \end{align}\tag{87}\] We can then rewrite (81 ) as \[\label{eq:l2hrd4a4} {\mathcal{L}}= \sum_{i=1}^{d} \left ({\boldsymbol{z}}_{1,i} {\boldsymbol{x}}_i^{(1)} + {\boldsymbol{z}}_{2,i}{\boldsymbol{x}}_i^{(2)} \right ) + \frac{1}{2} \sum_{i=1}^{d}\left ({\boldsymbol{x}}_i^{(1)}\right )^2{\boldsymbol{v}}_i + \frac{1}{2} \sum_{i=1}^{d}\left ({\boldsymbol{x}}_i^{(2)}\right )^2{\boldsymbol{v}}_i + \sum_{i=1}^{d}\nu_1{\boldsymbol{x}}_i^{(1)}{\boldsymbol{x}}_i^{(2)} -2\gamma -2\gamma_{0} c^2 -\nu q c^2.\tag{88}\] Set \[\begin{align} \label{eq:l2hrd4a5} K_{1,i} & = & 2\nu_i(-{\boldsymbol{z}}_{1,i}\nu_i + {\boldsymbol{z}}_{2,i}{\boldsymbol{v}}_i)({\boldsymbol{z}}_{1,i}{\boldsymbol{v}}_i - {\boldsymbol{z}}_{2,i}\nu_i) = 2 \nu_i(-z_{1,i}^2\nu_i{\boldsymbol{v}}_i + {\boldsymbol{z}}_{2,i}{\boldsymbol{z}}_{1,i}{\boldsymbol{v}}_i^2 + {\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i^2-{\boldsymbol{z}}_{2,i}^2\nu_i{\boldsymbol{v}}_i) \nonumber \\ K_{2,i} & = & {\boldsymbol{v}}_i(-{\boldsymbol{z}}_{1,i}\nu_i + {\boldsymbol{z}}_{2,i}{\boldsymbol{v}}_i)^2 = {\boldsymbol{v}}_i({\boldsymbol{z}}_{1,i}^2\nu_i^2 -2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i{\boldsymbol{v}}_i + {\boldsymbol{z}}_{2,i}^2{\boldsymbol{v}}_i^2) \nonumber \\ K_{3,i} & = & {\boldsymbol{v}}_i({\boldsymbol{z}}_{1,i}{\boldsymbol{v}}_i - {\boldsymbol{z}}_{2,i}\nu_i)^2 = {\boldsymbol{v}}_i({\boldsymbol{z}}_{1,i}^2{\boldsymbol{v}}_i^2 -2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i{\boldsymbol{v}}_i + {\boldsymbol{z}}_{2,i}^2\nu_i^2). \end{align}\tag{89}\] Then \[\label{eq:l2hrd4a6} {\mathcal{L}}= \sum_{i=1}^{d} \left ({\boldsymbol{z}}_{1,i} {\boldsymbol{x}}_i^{(1)} + {\boldsymbol{z}}_{2,i}{\boldsymbol{x}}_i^{(2)} \right ) + \frac{1}{2} \sum_{i=1}^{d} \frac{K_{1,i} + K_{2,i} + K_{3,i} }{(\nu_i^2-{\boldsymbol{v}}_i^2)^2} -2\gamma -2\gamma_{0} c^2 -\nu q c^2.\tag{90}\] From (89 ) we find \[\begin{align} \label{eq:l2hrd4a7} K_{1,i}+K_{2,i}+K_{3,i} & = ( -{\boldsymbol{z}}_{1,i}^2\nu_i^2{\boldsymbol{v}}_i + {\boldsymbol{v}}_i{\boldsymbol{z}}_{1,i}^2{\boldsymbol{v}}_i^2 ) +2(\nu_i^3{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i} -{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i{\boldsymbol{v}}_i^2 ) +(- 2{\boldsymbol{z}}_{2,i}^2\nu_i^2{\boldsymbol{v}}_i + {\boldsymbol{v}}_i{\boldsymbol{z}}_{2,i}^2\nu_i^2 +{\boldsymbol{v}}_i^3{\boldsymbol{z}}_{2,i}^2) \nonumber \\ & = {\boldsymbol{v}}_i{\boldsymbol{z}}_{1,i}^2( -\nu_i^2 + {\boldsymbol{v}}_i^2 ) +2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i( \nu_i^2 - {\boldsymbol{v}}_i^2 ) + {\boldsymbol{z}}_{2,i}^2{\boldsymbol{v}}_i(-\nu_i^2 + {\boldsymbol{v}}_i^2) \nonumber \\ & = ( -{\boldsymbol{v}}{\boldsymbol{z}}_{1,i}^2 + 2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i -{\boldsymbol{z}}_{2,i}^2{\boldsymbol{v}}_i )(\nu_i^2-{\boldsymbol{v}}_i^2). \end{align}\tag{91}\] Plugging this back in (90 ) gives \[\label{eq:l2hrd4a8} {\mathcal{L}}= \sum_{i=1}^{d} \left ({\boldsymbol{z}}_{1,i} {\boldsymbol{x}}_i^{(1)} + {\boldsymbol{z}}_{2,i}{\boldsymbol{x}}_i^{(2)} \right ) + \frac{1}{2} \sum_{i=1}^{d} \frac{ -{\boldsymbol{v}}_i{\boldsymbol{z}}_{1,i}^2 + 2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i -{\boldsymbol{z}}_{2,i}^2{\boldsymbol{v}}_i }{\nu_i^2-{\boldsymbol{v}}_i^2} -2\gamma -2\gamma_{0} c^2 -\nu q c^2.\tag{92}\] We also observe \[\begin{align} \label{eq:l2hrd4a9} {\boldsymbol{z}}_{1,i}{\boldsymbol{x}}_i^{(1)} + {\boldsymbol{z}}_{2,i}{\boldsymbol{x}}_i^{(2)} & = & {\boldsymbol{z}}_{1,i}\frac{ {\boldsymbol{z}}_{1,i}{\boldsymbol{v}}_i - {\boldsymbol{z}}_{2,i}\nu_i}{\nu_i^2-{\boldsymbol{v}}_i^2} + {\boldsymbol{z}}_{2,i}\frac{-{\boldsymbol{z}}_{1,i}\nu_i + {\boldsymbol{z}}_{2,i}{\boldsymbol{v}}_i}{\nu_i^2-{\boldsymbol{v}}_i^2} \nonumber \\ & = & \frac{ {\boldsymbol{z}}_{1,i}^2{\boldsymbol{v}}_i - 2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i + {\boldsymbol{z}}_{2,i}^2{\boldsymbol{v}}_i}{\nu_i^2-{\boldsymbol{v}}_i^2} . \end{align}\tag{93}\] A combination of (92 ) and (93 ) gives \[\label{eq:l2hrd4a10} \min_{{\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)}} {\mathcal{L}}= - \frac{1}{2} \sum_{i=1}^{d} \frac{ -{\boldsymbol{v}}_i{\boldsymbol{z}}_{1,i}^2 + 2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_i -{\boldsymbol{z}}_{2,i}^2{\boldsymbol{v}}_i }{\nu_i^2-{\boldsymbol{v}}_i^2} -2\gamma -2\gamma_{0} c^2 -\nu q c^2.\tag{94}\] After a change of variables \[\begin{align} \label{eq:l2hrd4a11} \nu & = &\nu_x\gamma_0 \nonumber \\ \nu_i & = &\nu_{x,i}\gamma_0 \nonumber \\ \nu_{x,i} & = &\nu_x{\boldsymbol{s}}_i^2 \nonumber \\ \gamma &=& \gamma_x\gamma_0 \nonumber \\ {\boldsymbol{v}}_{x,i} &=& 2(\gamma_x+{\boldsymbol{s}}_i^2) \nonumber \\ {\boldsymbol{v}}_i &= & 2(\gamma +\gamma_0{\boldsymbol{s}}_i^2)= 2\gamma_0(\gamma_x+{\boldsymbol{s}}_i^2) = {\boldsymbol{v}}_{x,i}\gamma_0, \end{align}\tag{95}\] (94 ) can be rewritten as \[\label{eq:l2hrd4a12} \min_{{\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)}} {\mathcal{L}}= - \frac{1}{4\gamma_0} \sum_{i=1}^{d} 2\frac{ -{\boldsymbol{v}}_{x,i}{\boldsymbol{z}}_{1,i}^2 + 2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_{x,i} -{\boldsymbol{z}}_{2,i}^2{\boldsymbol{v}}_{x,i} }{\nu_{x,i}^2-{\boldsymbol{v}}_{x,i}^2} - \gamma_0 ( 2\gamma_x +2 c^2 +\nu_x q c^2).\tag{96}\] Maximization ove \(\gamma_0\) gives \[\label{eq:l2hrd4a13} \max_{\gamma_0}\min_{{\boldsymbol{x}}^{(1)},{\boldsymbol{x}}^{(2)}} {\mathcal{L}}= - \sqrt{ \left (\sum_{i=1}^{d} 2\frac{ -{\boldsymbol{v}}_{x,i}{\boldsymbol{z}}_{1,i}^2 + 2{\boldsymbol{z}}_{1,i}{\boldsymbol{z}}_{2,i}\nu_{x,i} -{\boldsymbol{z}}_{2,i}^2{\boldsymbol{v}}_{x,i} }{\nu_{x,i}^2-{\boldsymbol{v}}_{x,i}^2} \right )( 2\gamma_x +2 c^2 +\nu_x q c^2) }.\tag{97}\] Combining (82 ), (95 ), and (97 ) and relying on the law of large numbers and concentrations, we obtain \[\begin{align} \label{eq:l2hrd14} \lim_{d\rightarrow \infty }\frac{1}{\sqrt{d}} {\mathbb{E}}L^{(2)}(c;t) & \leq & \lim_{d\rightarrow \infty } \min_{\gamma_x,\nu_x} \sqrt{ \left (\frac{1}{d} \sum_{i=1}^{d} 2{\boldsymbol{s}}_i^2\frac{ -2{\boldsymbol{v}}_{x,i} + 2a\nu_{x,i} }{\nu_{x,i}^2-{\boldsymbol{v}}_{x,i}^2} \right )( 2\gamma_x +2 c^2 +\nu_x q c^2) } \nonumber \\ & = & \lim_{d\rightarrow \infty } \min_{\gamma_x,\nu_x} \sqrt{ \left (\frac{1}{d} \sum_{i=1}^{d} 4{\boldsymbol{s}}_i^2\frac{ {\boldsymbol{v}}_{x,i} - a\nu_{x,i} }{{\boldsymbol{v}}_{x,i}^2 - \nu_{x,i}^2} \right )( 2\gamma_x +2 c^2 +\nu_x q c^2) } \nonumber \\ & = & \lim_{d\rightarrow \infty } \min_{\gamma_x,\nu_x} \sqrt{ \left (\frac{1}{d} \sum_{i=1}^{d} 4{\boldsymbol{s}}_i^2\frac{ 2(\gamma_x+{\boldsymbol{s}}_i^2) - a\nu_x{\boldsymbol{s}}_i^2 }{4(\gamma_x+{\boldsymbol{s}}_i^2)^2 - \nu_x^2{\boldsymbol{s}}_i^4} \right )( 2\gamma_x +2 c^2 +\nu_x q c^2) }. \end{align}\tag{98}\]

The following theorem summarizes the above characterization of \(\lim_{d\rightarrow \infty }\frac{1}{\sqrt{d}} {\mathbb{E}}L^{(2)}(c;t)\).

Theorem 7 provides all needed ingredients so that one can numerically check whether (?? ) holds. As noted in Remark 4, when the true covariance \(\Sigma\) is random and characterized via its spectral density, all empirical averages become integrals. One then numerically evaluates \(\bar{L}^{(2)}\) and \(\bar{L}\) and if (?? ) holds, the proof is completed up to the level of numerical precision. Since numerical precision can be arbitrary one can effectively approach equality in (?? ) arbitrarily closely. We tested quite a few ensembles and always obtained that (?? ) holds.

It is even more interesting to apply the same methodology to deterministic setup presented in the previous section. In that case one can not numerically rely on infinite sizes. Nonetheless, one can expect the very same behavior for sufficiently large \(d\). To see whether this indeed happens, we conducted a set of numerical evaluations and the obtained results are shown in Figure 1. We selected \(d=3000\) and \({\boldsymbol{s}}=linspace_d[0.5,1]\) (\(linspace_d[0.5,1]\) is an increasing vector of \(d\) components equally spaced between \(0.5\) and \(1\); in other words, \({\boldsymbol{s}}_i-{\boldsymbol{s}}_{i-1}={\boldsymbol{s}}_j-{\boldsymbol{s}}_{j-1}\) for any \(i\) and \(j\); both endpoints of the interval are included as well, i.e., \({\boldsymbol{s}}_1=0.5\) and \({\boldsymbol{s}}_d=1\)). While in theory \(d\rightarrow \infty\), we found that \(d\) on the order of a few thousands already behaves very well. As Figure 1 indicates, there are no local optima and the condition (?? ) from Theorem 7 is indeed met (none of the curves reaches the maximum except at \(q=1\)).

We complement the above numerical findings with the following theorem which analytically confirms that condition (?? ) is met.

3.5 Matching \(\delta\) and \({\mathbb{E}}\lambda_n(\hat{\Sigma}-\Sigma)\)↩︎

We first recall/summarize all key aspects of the above analysis. From (13 ), we have \[\begin{align} \label{eq:dellam1} \lambda_n\left (\hat{\Sigma} -\Sigma \right ) & = & \max_{c} \left (\left (\xi(c)\right )^2 - c^2 \right ). \end{align}\tag{110}\] Theorems 1 and 4 (equations (?? ) and (?? )) give \[\begin{align} \label{eq:dellam2} c + \frac{1}{\sqrt{2n}c} {\mathbb{E}}B(c) \leq {\mathbb{E}}\xi(c) \leq c + \frac{1}{\sqrt{n}} {\mathbb{E}}L(c). \end{align}\tag{111}\] Also, Theorems 7 and 8 give \[\begin{align} \label{eq:dellam3} \lim_{d\rightarrow \infty } \frac{1}{\sqrt{2d}c}{\mathbb{E}}B(c) = \lim_{d\rightarrow \infty } \frac{1}{\sqrt{d}} {\mathbb{E}}L(c), \end{align}\tag{112}\] and \[\begin{align} \label{eq:dellam4} \lim_{d\rightarrow \infty } {\mathbb{E}}\xi(c) = c + \frac{1}{\sqrt{\alpha}}\lim_{d\rightarrow \infty } \frac{1}{\sqrt{d}} {\mathbb{E}}L(c). \end{align}\tag{113}\] Combining (110 ) and (113 ), we further find \[\begin{align} \label{eq:dellam5} \lim_{d\rightarrow \infty } {\mathbb{E}}\lambda_n\left (\hat{\Sigma} -\Sigma \right ) & = & \max_c \left (\frac{2c}{\sqrt{\alpha}}\lim_{d\rightarrow \infty } \frac{1}{\sqrt{d}} {\mathbb{E}}L(c) + \frac{1}{\alpha}\left (\lim_{d\rightarrow \infty } \frac{1}{\sqrt{d}} {\mathbb{E}}L(c) \right )^2 \right ). \end{align}\tag{114}\] The above is already a very convenient characterization of sample covariance estimation error. It basically relates to the maximal eigenvalue of a residual error matrix \(\hat{\Sigma}-\Sigma\). One can go a step further and show that this also matches the corresponding characterization of the spectral norm as well.

In particular, from (8 ), we have \[\label{eq:delam7} \|\hat{\Sigma} - \Sigma\|_2 = \max(\lambda_n(\hat{\Sigma} - \Sigma),|\lambda_1(\hat{\Sigma} - \Sigma)|).\tag{115}\] The discussion from previous sections ensures that the righthand side of (113 ), i.e., \(\lim_{d\rightarrow \infty } {\mathbb{E}}\lambda_n ( \hat{\Sigma} -\Sigma )\), is positive. If \(\lim_{d\rightarrow \infty } {\mathbb{E}}\lambda_1(\hat{\Sigma} - \Sigma)\geq 0\) one then automatically has \[\label{eq:delam9} \lim_{d\rightarrow \infty } {\mathbb{E}}\|\hat{\Sigma} - \Sigma\|_2 = \lim_{d\rightarrow \infty } {\mathbb{E}}\lambda_n(\hat{\Sigma} - \Sigma) .\tag{116}\] Therefore, the only interesting remaining scenario to consider is \(\lim_{d\rightarrow \infty } {\mathbb{E}}\lambda_1 (\hat{\Sigma} - \Sigma) < 0\). To that end, we observe that when \(\lambda_1 (\hat{\Sigma} - \Sigma) < 0\) then \[\begin{align} \label{eq:dellam10} \left |\lambda_1 \left (\hat{\Sigma} - \Sigma\right )\right | = \lambda_n \left (\Sigma -\hat{\Sigma}\right ) \end{align}\tag{117}\] Following (9 )–(11 ) \[\begin{align} \label{eq:dellam11} \lambda_n \left (\Sigma -\hat{\Sigma}\right )& = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d} \left ({\boldsymbol{x}}^T SS^T {\boldsymbol{x}}-\frac{1}{n}{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}\right ). \end{align}\tag{118}\] After setting \[\begin{align} \label{eq:dellam12} \xi^-(c) & = & \frac{1}{\sqrt{n}} \min_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c } \sqrt{{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}} = \frac{1}{\sqrt{n}} \min_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c }\max_{{\boldsymbol{y}}\in{\mathbb{S}}^n } {\boldsymbol{y}}AS {\boldsymbol{x}}, \end{align}\tag{119}\] we have in place of (118 ) \[\begin{align} \label{eq:dellam13} \lambda_n \left (\Sigma -\hat{\Sigma}\right )& = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d} \left ({\boldsymbol{x}}^T SS^T {\boldsymbol{x}}- \frac{1}{n}{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}\right ) \nonumber \\ & = & \max_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c,c} \left ({\boldsymbol{x}}^T SS^T {\boldsymbol{x}}- \frac{1}{n}{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}\right ) \nonumber \\ & = & \max_{c} \left (c^2 - \min_{{\boldsymbol{x}}\in{\mathbb{S}}^d,\|S{\boldsymbol{x}}\|_2=c} \frac{1}{n}{\boldsymbol{x}}^T S^TA^TAS {\boldsymbol{x}}\right ) \nonumber \\ & = & \max_{c} \left (c^2 - \left (\xi^-(c)\right )^2 \right ). \end{align}\tag{120}\] To characterize \(\xi^-(c)\) and consequently \(\lambda_n\left (\Sigma - \hat{\Sigma} \right )\), we again utilize RDT. The following theorem summarizes the obtained results.

Finally, we are now in position to complete the story and establish spectral radius analogue to Theorems 3 and 7.

3.6 Practical numerical utilization↩︎

The results from the previous sections provide characterization of \(\lambda_n\left (\hat{\Sigma} -\Sigma \right )\) and \(\delta(\alpha)\). In Figure 2 we compare these theoretical predictions with simulated results. As in Figure 1, for different dimensions \(d\), we selected \({\boldsymbol{s}}=linspace_d[0.5,1]\) as a vector of \(d\) equally spaced components from \([0.5,1]\) interval. While the theoretical framework requires \(d\rightarrow \infty\), Figure 2 shows that an excellent agreement between the theory and simulations is achieved already for \(d\) on the order of thousand.

The conducted analysis allows to obtain very precise estimation error characterizations. Consequently, it enables to accurately predict concrete effects of the increased number of samples. We show this in Figure 3 where \(\alpha\) is successively increased. As can be seen, the simulated results again closely follow the theoretical predictions, although the utilized dimensions are relatively small compared to the infinitely-dimensional setup required by the theoretical analysis.

4 Conclusion↩︎

We studied the sample covariance error of centered Gaussians. To move beyond scaling characterizations and determine the precise limiting value of the error’s spectral norm, we have developed a generic framework based on Random Duality Theory (RDT). The framework consists of three key components: (1) Deriving closed-form, explicit RDT-based upper bounds; (2) Introducing a novel bilinear-quadratic RDT lower-bounding mechanism; and (3) Combining this lower-bounding mechanism with a 2-replica systems bounding strategy to demonstrate that the upper bounds can be matched in large-dimensional contexts. Our theoretical predictions show excellent agreement with the results obtained through numerical evaluations and simulations.

The developed framework is highly generic and offers significant potential for further extensions, including studying nearly all associated problem variants considered in recent literature and beyond. While the conceptual foundations remain the same, the specific underlying technical considerations are problem-dependent and will be discussed in future work.

References↩︎