The Quadratic Optimization Bias of Large Covariance Matrices

To motivate the study of the interplay between optimization and model estimation error, we consider a quadratic function in $p$ variables. Let,

for an inner product $\langle\hskip 1.0pt\cdot,\cdot\hskip 1.0pt\rangle$, constants $c_{0},c_{1}\in\bbR$ and a vector $\zeta\in\bbR^{p}$. The $p\times p$ matrix $\Sigma$ is assumed to be symmetric and positive definite. The maximization of $(\ref{quad})$ is encountered in many classical contexts within statistics and probability including least-squares regression, maximum a posteriori estimation, saddle-point approximations, and Legendre-Fenchel transforms in moderate/large deviations theory. Some related and highly influential applications include Markowitz’s portfolio construction in finance (Markowitz, 1952), Capon beamforming in signal processing (Capon, 1969) and optimal fingerprinting in climate science (Hegerl et al., 1996). In optimization theory, quadratic functions form a key ingredient for more general (black-box) minimization techniques such as trust-region methods (e.g., Maggiar et al. (2018)).¹¹1In this setting the covariance matrix corresponds to an estimated Hessian matrix. Since any number of linear equality constraints may be put into the unconstrained Lagrangian form in $(\ref{quad})$, our setting is more general than it first appears. Moreover, the maximization of $(\ref{quad})$ is the starting point for numerous applications of quadratic programming where nonlinear constraints are often added.²²2To give one example, interpreting $\Sigma$ as a graph adjacency matrix and adding simple bound constants to $(\ref{quad})$ leads to approximations of graph properties such as the maximum independent set (Hager and Hungerford, 2015). While $\Sigma$ is no longer interpreted as a spiked covariance matrix in a graph theory setting, its mathematical properties are similar owing to the celebrated Cheeger’s inequality.

The maximizer of $Q(\hskip 1.0pt\cdot\hskip 1.0pt)$ is given by $c_{1}\hskip 1.0pt\Sigma^{-1}\zeta$ which attains the objective value

but in practice, the maximization of $Q(\hskip 1.0pt\cdot\hskip 1.0pt)$ is performed with an estimate $\bm{\hat{}}{\Sigma}$ replacing the unknown $\Sigma$. This “plug-in” step is well known to yield a perplexing problem (see Section 1.3). In essence, the optimizer chases the errors in $\bm{\hat{}}{\Sigma}$ to produce a systematic bias in the computed maximum. This bias is then amplified by a higher dimension.

Consider a high-dimensional limit $p\uparrow\infty$ and a sequence of symmetric positive definite $\Sigma=\Sigma_{p\times p}$ with a fixed number $q$ of spiked eigenvalues diverging in $p$ and all remaining eigenvalues bounded in $(0,\infty)$. Let $\hat{x}$ be the maximizer of $\hat{Q}(\hskip 1.0pt\cdot\hskip 1.0pt)$, defined by replacing $\Sigma$ in $(\ref{quad})$ by estimates $\bm{\hat{}}{\Sigma}=\bm{\hat{}}{\Sigma}_{p\times p}$ with the same eigenvalue properties. The estimated objective is $\hat{Q}(\hat{x})$, but a more relevant quantity is the realized objective,

where $\hat{\mu}^{2}_{p}=\langle\zeta,\bm{\hat{}}{\Sigma}^{-1}\zeta\rangle$ and $\hat{D}_{p}$ is a discrepancy (relative to $(\ref{conv})$) that can grow rapidly as the dimension increases. Precluding edge cases where $\hat{\mu}^{2}_{p}\hskip 1.0pt/\langle\zeta,\zeta\rangle$ or $\langle\zeta,\zeta\rangle$ vanish, unless $\bm{\hat{}}{\Sigma}$ is fine-tuned in a calculated way, the following puzzling behavior ensues.

1.  

   The discrepancy $\hat{D}_{p}$ tends to $-\infty$ as $p\uparrow\infty$ and consequently, the realized maximum $Q(\hat{x})$ tends to $-\infty$ while the true maximum $(\ref{conv})$ tends to $+\infty$.

The asymptotic behavior above is fully determined by a certain $\bbR^{q}$-valued function $\scrE_{p}(\hskip 1.0pt\cdot\hskip 1.0pt)$ which we derive and call the quadratic optimization bias. The way in which this bias depends on the entries of $\bm{\hat{}}{\Sigma}$ characterizes the sought after interplay between the optimizer and the error in the estimated covariance model. Mitigating the discrepancy between the realized and true quadratic optima $(\ref{conv})$ and $(\ref{Qhx})$ reduces to the problem of approximating the roots of $\scrE_{p}(\hskip 1.0pt\cdot\hskip 1.0pt)$. We remark that by parametrizing the constants $c_{0}$ and $c_{1}$ in $p$, one can arrive at an alternative limits for $(\ref{conv})$ and $(\ref{Qhx})$, but practical scalings preserve the large disparity between the true and realized objective values. We examine (in Section 2.2) the scaling $c_{1}=1/p$ in particular, due to its applicability to portfolio theory, robust beamforming and optimal fingerprinting.

The illustration above reflects that, in statistical settings, solutions to estimated quadratic optimization problems exhibit very poor properties out-of-sample. Section 2 answers the question of which components of $\Sigma$ must be estimated accurately to reduce the discrepancy $\hat{D}_{p}$ in $(\ref{Qhx})$. The size of $|\hat{D}_{p}|$ is amplified by the growth rate $r_{p}$ of the $q$ spiked eigenvalues, but is fully determined by the precision of the estimate $\scrH$ of the associated $p\times q$ matrix of eigenvectors $\scrB$ of $\Sigma$. In particular, $\hat{D}_{p}=-|\scrE_{p}(\scrH)|^{2}\hskip 1.0ptO(r_{p})$ where $\scrE_{p}(\scrH)$ is given by,

for the Euclidean length $|\cdot|\hskip 1.0pt$. Theorem 1 gives sharp asymptotics for $\hat{D}_{p}$ in $\scrE_{p}(\scrH)$ and the other estimates/parameters. Remarkably, the accuracy of the estimates of eigenvalues of $\Sigma$ is secondary relative the ensuring that $\scrH$ is such that $\scrE_{p}(\scrH)$ is small for large $p$. This is noteworthy in view of the large literature on bias correction of sample eigenvalues (or “eigenvalue shrinkage”: see Ollila, Palomar and Pascal (2020), Ledoit and Wolf (2021), Ledoit and Wolf (2022) and Donoho, Gavish and Romanov (2023) for a sampling of recent work). Instead, for the discrepancy $\hat{D}_{p}$, the estimation of the eigenvectors of the spikes is what matters most. We remark that while $\scrH=\scrB$ forms a root of the map $\scrE_{p}\hskip 0.5pt:\hskip 0.5pt\bbR^{p\times q}\to\bbR^{q}$ (i.e., $\scrE_{p}(\scrB)=0_{q}$), it is not the only root. We refer to $\scrE_{p}(\hskip 1.0pt\cdot\hskip 1.0pt)$ as the quadratic optimization bias (function) which was first identified in Goldberg, Papanicolaou and Shkolnik (2022) in the context of portfolio theory and for the special covariance $\Sigma$ with a single spike $(q=1)$ and identical remaining eigenvalues.³³3We state a more general definition in Section 2, but $\scrH$ must have orthonormal columns in $(\ref{obf})$.

Section 3 considers a $p\times p$ sample covariance matrix $S$ and its spectral decomposition $S=\scrH\calS_{p}^{2}\scrH^{\top}+G$, for a diagonal $q\times q$ matrix of eigenvalues $\calS^{2}_{p}$, the associated $p\times q$ matrix $\scrH$ of eigenvectors ($\scrH^{\top}\scrH=I_{q}$) and a residual $G$. It is assumed that $r_{p}$ is $O(p)$ and that the sequence $S=S_{p\times p}$ is based on a fixed number of observations of a high dimensional vector. Our Theorem 3 proves that $\scrE_{p}(\scrH)$ is almost surely bounded away from zero (in $\bbR^{q}$) eventually in $p$. This has material implications for the use of principal component analysis for problems motivated by Section 1.1.

Section 5 develops the following correction to the sample eigenvectors $\scrH$. For the $q\times q$ diagonal matrix ${{\Psi}}$ satisfying ${{\Psi}}^{2}=I_{q}-\hskip 1.0pt\textup{\text{tr}}(G)\hskip 1.0pt\calS_{p}^{-2}/n_{q}$ for $n_{q}\geq 1$, the difference between the number of nonzero sample eigenvalues and $q$, we compute

Theorem 9 proves the $p\times q$ matrix of left singular vectors of $(\ref{Hcorr})$, denoted $\scrH_{\sharp}$, has

almost surely. The matrix $\scrH_{\sharp}$ constitutes a set of corrected principal component loadings and is the basis of our covariance estimator $\bm{\hat{}}{\Sigma}_{\sharp}$. This matrix, owing to $(\ref{EHslim})$, yields an improved plug-in estimator $x_{\sharp}=c_{1}\bm{\hat{}}{\Sigma}_{\sharp}^{-1}\zeta$ for the maximizer of $(\ref{quad})$. Thus, our work also has implications for the estimation of the precision matrix $\Sigma^{-1}$. Theorem 9 also proves that the columns of $\scrH_{\sharp}$ have a larger projection (than $\scrH$) onto the column space of $\scrB$. Recent literature has remarked on the difficulty (or even impossibility) of correcting such bias in eigenvectors (e.g., Ledoit and Wolf (2017), Wang and Fan (2017) and Jung (2022)). That projection is strictly better when $z$ in $(\ref{obf})$ has $|\scrB^{\top}z|$ bounded away from zero, i.e., captures information about that subspace. But $(\ref{EHslim})$ holds regardless, highlighting that the choice of the “loss” function (in our case $(\ref{Qhx})$) matters.⁴⁴4See also Donoho, Gavish and Johnstone (2018) for another illustration of this phenomenon.

In Section 4, we prove an impossibility theorem (Theorem 8) that shows that without very strong assumptions one cannot obtain an estimator of $\scrB^{\top}\scrH$ asymptotically in the dimension if $q>1$. This has negative implications for obtaining estimates of $\scrE_{p}(\scrH)$ in $(\ref{obf})$ where $\scrB^{\top}\scrH$ is one of the unknowns. The latter contains all $q^{2}$ inner products between the sample and population eigenvectors, and its estimation from the observed data is an interesting theoretical problem in its own right. Our negative result adds to the literature on high dimension and low sample size (hdlss) asymptotics, as inspired by Hall, Marron and Neeman (2005) and Ahn et al. (2007).⁵⁵5Aoshima et al. (2018) survey much of the literature since. We remark that the hdlss regime is highly relevant for real-world data as a small sample size is often imposed by experimental constraints, or by the lack of long-range stationarity of time series. The content of Theorem 8 also points to a key feature that distinguishes our work from Goldberg, Papanicolaou and Shkolnik (2022)) who fix $q=1$. Another aspect making our setting substantially more challenging is that we find roots of a multivariate function $\scrE_{p}(\hskip 1.0pt\cdot\hskip 1.0pt)$ (which is univariate when $q=1$).

In terms of applications, our results generalize those of Goldberg, Papanicolaou and Shkolnik (2022) to covariance models that hold wide acceptance in the empirical literature on financial asset return (i.e., the Arbitrage Pricing Theory of Ross (1976), Huberman (1982), Chamberlain and Rothschild (1983) and others). Section 6 investigates the problem of minimum variance investing with numerical simulation, and demonstrates that the estimator $\scrH_{\sharp}$ results in vanishing asymptotic portfolio risk and a bounded discrepancy $\hat{D}_{p}$ (see Figure 1). Appendix E summarizes other applications including signal-noise processing and climate science as related to Section 1.1.

Our findings in Section 1.1 form a starting point for important extensions and applications. Extending the estimator in $(\ref{Hcorr})$ to general quadratic programming with inequality constraints would greatly expand its scope. In terms of covariance models, we require spikes that diverge linearly with the dimension, which excludes several alternative frameworks in the literature.⁶⁶6This includes the Johnstone spike model, in which all eigenvalues remain bounded as the dimension grows, and its extensions (e.g., Johnstone (2001), Paul (2007), Johnstone and Lu (2009) and Bai and Yao (2012)). Futher generalizations include slowly growing spiked eigenvalue models as in De Mol, Giannone and Reichlin (2008), Onatski (2012), Shen et al. (2016) and Bai and Ng (2023). Likewise, the asymptotics of the data matrix aspect ratio $p/n$ differs across applications. We also do not address the important setting in which the number of spikes $q$ is misspecified.⁷⁷7There is a large literature on the estimation of the number of spikes/factors/principal components. Most relevant to our setup (high dimension and low sample size) is Jung, Lee and Ahn (2018). Finally, the established convergence in $(\ref{EHslim})$ leaves the question of rates unanswered. This is particularly important for problems requiring the discrepancy $\hat{D}_{p}$ to not grow too quickly. We offer no theoretical treatment of convergence rates but our numerical results suggest this quantity remains bounded as $p$ grows (c.f., Figure 1).

The work we build on directly was initiated in Goldberg, Papanicolaou and Shkolnik (2022). We refer to their proposal as the GPS estimator and derive it in Section 5.1. Important extensions are developed in Gurdogan and Kercheval (2022) and Goldberg, Gurdogan and Kercheval (2023). The GPS estimator was shown to be mathematically equivalent to a James-Stein estimation of the leading eigenvector of a sample covariance matrix in Shkolnik (2022). These results share much in common with the ideas found in Casella and Hwang (1982). For a survey of the above literature, focusing on connections to the James-Stein estimator, see Goldberg and Kercheval (2023). The GPS estimator is explained in terms of regularization in Lee and Shkolnik (2024a), and Lee and Shkolnik (2024b) derive central limit theorems for this estimator as relevant for the convergence of the discrepancy $\hat{D}_{p}$. Some numerical exploration of the case of more than one spike is found in Goldberg et al. (2020).

The spiked covariance models we consider, and the application of pca for their estimation, are rooted in the literature on approximate factor models and “asymptotic principal components” originating with Chamberlain and Rothschild (1983) and Connor and Korajczyk (1986). Recent work in this direction is well represented by Bai and Ng (2008), Fan, Liao and Mincheva (2013), Bai and Ng (2023) and Fan, Masini and Medeiros (2023). In this literature, the work that most closely resembles ours, by focusing on improved estimation of sample eigenvectors, is Fan, Liao and Wang (2016), Fan and Zhong (2018) and Lettau and Pelger (2020). Fan, Liao and Wang (2016) project the data onto a space generated by some externally observed covariates, improving the resulting sample eigenvectors when the covariates have sufficient explanatory power. Fan and Zhong (2018) apply a linear transformation to the sample eigenvectors in an approach that is most closely related to formula $(\ref{Hcorr})$. We also apply a linear transformation, but the eigenspace is first augmented by the vector $\zeta$ in $(\ref{quad})$.⁸⁸8We remark that with a single spike/factor (i.e., $q=1$), a linear transformation of the eigenvector(s) adjustment only the eigenvalue, not the eigenvector itself due to its unit length normalization. Further differences with Fan, Liao and Wang (2016) arise in the estimation of the optimal linear transformation. Lettau and Pelger (2020) extract principal components from a rank-one updated sample covariance matrix. This update is based on insight from asset pricing theory and it is unclear how the resulting sample eigenvectors are related to formula $(\ref{Hcorr})$. The same applies to the very closely related literature on sample covariance matrix shrinkage (e.g., Ledoit and Wolf (2004a), Fisher and Sun (2011), Lancewicki and Aladjem (2014) and Wang and Zhang (2024)).⁹⁹9This takes the form $\bm{\hat{}}{\Sigma}=a\hskip 1.0ptS+(1-a)F$ for some $a\in[0,1]$ and matrix $F$. Targets $F\neq I$ adjust eigenvectors but in ways that may be difficult to quantify via closed-form expressions (c.f. $(\ref{Hcorr})$).

The vast majority of the literature on approximate factor models and covariance estimation assumes the data matrix aspect ratio $p/n$ tends to a finite constant asymptotically.¹⁰¹⁰10This may be due to the outsized influence of random matrix theory (e.g., Marchenko and Pastur (1967)). Another reason may be the consistency of the sample eigenvectors that can be achieved in this regime (see Yata and Aoshima (2009), Shen, Shen and Marron (2016) and Wang and Fan (2017). In contrast, our analysis of a finite sample in the high dimensional limit draws on the work on pca in Jung and Marron (2009), Jung, Sen and Marron (2012) and Shen et al. (2016) and others. In the latter, the hdlss asymptotics for the matrix $\scrB^{\top}\scrH$, appearing in $(\ref{obf})$, have already been worked out (but see Section 4 for our impossibility theorem). Our main focus is on correcting the biases that the asymptotics of $\scrB^{\top}\scrH$ reveal. For approaches to correcting the finite sample bias in eigenvalues and principal component scores, see Yata and Aoshima (2012), Yata and Aoshima (2013), Jung (2022) and our Remark 7. Shen, Shen and Marron (2013) apply regularization in the presence of sparsity in the population eigenvectors to correct finite sample bias in the principal components. It is unclear how their estimators are related to the update in $(\ref{Hcorr})$, but we do not impose such sparsity assumptions.

Several other strands of the pca literature are relevant as their aims coincide with improved sample eigenvector estimation. In one direction is the literature on sparse and low-rank matrix decompositions (e.g. Chandrasekaran, Parrilo and Willsky (2012), Saunderson et al. (2012), Bai and Ng (2019), Farnè and Montanari (2024) and Li and Shkolnik (2024)). These convex relaxations aim to find more accurate low-dimensional representations of the data and are sometime referred to as forms of robust pca (Candès et al., 2011). In a related direction is the recent work on robust pca for heteroskedastic noise (e.g., Cai et al. (2021), Zhang, Cai and Wu (2022), Yan, Chen and Fan (2021), Agterberg, Lubberts and Priebe (2022) and Zhou and Chen (2023)). These efforts provide ($p,n$ finite) bounds on generalized angles between the true and the estimated subspaces and complement our asymptotic pca results in Sections 3 & 5. Perturbations of eigenvectors have also been recently revisited in Fan, Wang and Zhong (2018), Abbe, Fan and Wang (2022) and Li et al. (2022). The latter use these bounds to construct estimators that “de-bias” linear forms such as $\scrB^{\top}z$ appearing in $(\ref{obf})$. These results can likely supply alternative proofs to ours (or even convergence rates), but our focus is on limit theorems only.

Lastly, we emphasize the area of mean-variance portfolio optimization. As the literature on this topic is quite vast, we mention only a few strands related to Section 1.1. Examples of early influential work in this direction include Michaud (1989) and Best and Grauer (1991). For numerical simulations that illustrate the impact on practically motivated models and metrics see Bianchi, Goldberg and Rosenberg (2017). A random matrix theory perspective on the behavior of objectives related to $(\ref{Qhx})$ may be found in Pafka and Kondor (2003), Bai, Liu and Wong (2009), El Karoui (2010), El Karoui (2013), Bun, Bouchaud and Potters (2017) and Bodnar, Okhrin and Parolya (2022). Highly relevant recent work in econometrics using latent factor models includes Ding, Li and Zheng (2021) who consider a portfolio risk measure closely tied to $(\ref{Qhx})$. Bayesian approaches to mean-variance optimization include Lai et al. (2011) and Bauder et al. (2021). These estimators are closely related to Ledoit-Wolf shrinkage (Ledoit and Wolf (2003) and Ledoit and Wolf (2004b)) which itself has undergone numerous improvements (e.g., Ledoit and Wolf (2018) and Ledoit and Wolf (2020a)). In tandem, shrinkage methods have been known to impart effects akin to extra constraints in the portfolio optimization as early as Jagannathan and Ma (2003). An insightful example of such robust portfolio optimization that relates $(\ref{Qhx})$ to the convergence of the covariance matrix estimator is developed in Fan, Zhang and Yu (2012). More advanced robust portfolio optimizations have also been proposed (e.g., Boyd et al. (2024)). Alternatively, constraints are often applied in the covariance matrix estimation process as an optimization in itself. For example, Won et al. (2013) apply a condition number constraint that leads to non-linear adjustments of sample eigenvalues (c.f., Ledoit and Wolf (2020b)), but leaves the sample eigenvectors unchanged. Bongiorno and Challet (2023) document the difficulty with relying solely on eigenvalue correction, especially for small sample sizes. Cai et al. (2020) apply sparsity constraints (on the precision matrix) and analyze optimality properties related to $(\ref{Qhx})$. We emphasize that the impact of such constraints on eigenvectors is difficult (or impossible) to quantify, in contrast to formula $(\ref{Hcorr})$.¹¹¹¹11It should be noted that another interesting approach to mean-variance portfolio optimization concerns the direct shrinkage of the portfolio weights (i.e., akin to shrinkage $\hat{x}$ in $(\ref{Qhx})$, e.g., Bodnar, Parolya and Schmid (2018), Bodnar, Okhrin and Parolya (2022) Bodnar, Parolya and Thorsén (2023)).

Let $\textsc{col}{\hskip 0.5pt(A)}$ denote the column span of the matrix $A$ and let $\zeta_{A}$ denote the orthogonal projection of the vector $\zeta$ on $\textsc{col}{\hskip 0.5pt(A)}$, e.g.,

where $A^{\dagger}=(A^{\top}A)^{-1}A^{\top}$, the Moore-Penrose inverse of a full column rank matrix $A$. We use $I$ to denote an identity matrix and $I_{q}$ when highlighting its dimensions, $q\times q$.

Write $\langle u,v\rangle$ for the scalar product of $u,v\in\bbR^{m}$, let $|u|=\sqrt{\langle u,u\rangle}$ and $|A|$ be the induced (spectral) norm of a matrix $A$. We denote by $\nu_{m\times q}(\hskip 1.0pt\cdot\hskip 1.0pt)$ a function that given a $m\times\ell$ matrix $A$, uniquely selects (see Appendix D) an enumeration of its singular values $|A|=\Lambda_{11}\geq\cdots\geq\Lambda_{\min(\ell,m)}$ and outputs a $m\times q$ matrix $\nu_{m\times q}(A)$ of left singular vectors with the values $\Lambda_{11},\dots,\Lambda_{qq}$ in columns $1,\dots,q\leq\min(m,\ell)$. That is,

where $\Lambda$ is the $q\times q$ diagonal with entries $\Lambda_{11},\dots,\Lambda_{qq}$, and $\scrR$ is the $\ell\times q$ matrix of right singular vectors of $A$. The $m\times q$ matrix $\nu_{m\times q}(A)$ also corresponds to some unique choice of eigenvectors of $AA^{\top}$ with $q$ largest eigenvalues $\Lambda_{11}^{2}\geq\cdots\geq\Lambda^{2}_{qq}$.

We take $0_{q}=(0,\dots,0)\in\bbR^{q}$ and $1_{q}=(1,\dots,1)\in\bbR^{q}$. Lastly, $\varliminf_{p\uparrow\infty}$ and $\varlimsup_{p\uparrow\infty}$ denote the limit inferior and superior as $p\uparrow\infty$, and $A=A_{p\times m}$, a sequence of matrices with dimensions $p\times m$ when at least one of $p$ or $m$ grows to infinity.

Returning to the optimization setting of Section 1.1, for constants $c_{0},c_{1}\in\bbR$ and $\zeta\in\bbR^{p}$, we consider

which attains $\max_{x\in\bbR^{p}}\hat{Q}(x)=\hat{Q}(\hat{x})=c_{0}+\frac{c_{1}^{2}\hat{\mu}_{p}^{2}}{2}$ at the maximizer $\hat{x}\in\bbR^{p}$ analogously to $(\ref{conv})$ but with $\hat{\mu}_{p}^{2}=\langle\zeta,\bm{\hat{}}{\Sigma}^{-1}\zeta\rangle$. Because $\hat{Q}(\hskip 1.0pt\cdot\hskip 1.0pt)$ is not the true objective function $Q(\hskip 1.0pt\cdot\hskip 1.0pt)$ in $(\ref{quad})$, we are interested in the realized objective $Q(\hat{x})$. Now,

which identifies the discrepancy $\hat{D}_{p}$ in $(\ref{Qhx})$ relative to both $\hat{Q}(\hat{x})$ and $(\ref{conv})$.

To avoid division by zero in $(\ref{optbias})$, we prevent $\zeta\in\bbR^{p}$ from vanishing and residing entirely in $\textsc{col}{\hskip 0.5pt(H)}$ asymptotically (i.e., $|1-z_{H}|^{2}=1-|\zeta_{H}|^{2}/|\zeta|^{2}$).¹³¹³13This edge case must be treated separately from our analysis and we do not pursue it. The entries of $\zeta$ may be viewed as the first $p$ entries of an infinite sequence or as rows of a triangular array. We further assume the estimate $H$ has properties consistent $B$ in view of Assumption 1 (b).

We address the asymptotics of the discrepancy $\hat{D}_{p}=2-\langle\hat{x},\Sigma\hat{x}\rangle/(c_{1}\hat{\mu}_{p})^{2}$ in $(\ref{QhxD})$, letting $BB^{\top}=\scrB\Lambda_{p}^{2}\scrB$ as above, with $\scrB=\nu_{p\times q}(B)$ the canonical choice.

Theorem 1 reveals that $\hat{D}_{p}$ diverges to $-\infty$ unless we find roots of $\scrE_{p}(\hskip 1.0pt\cdot\hskip 1.0pt)$, perhaps asymptotically. Note that $\scrE_{p}(B)=0$, but other roots exists (see Section 5).

Lemma 2 pinpoints what constitutes a poor “plug-in” covariance estimator $\bm{\hat{}}{\Sigma}$. For example, the column lengths of $H$ have no effect on the quadratic optimization bias $\scrE_{p}(H)$. For the eigenvalue decomposition $HH^{\top}=\scrH\calS_{p}^{2}\scrH^{\top}$ (with $K=\calS_{p}^{-1}$ in Lemma 2), we see that $\scrE_{p}(H)=\scrE_{p}(\scrH)$. Thus, to fine-tune $\bm{\hat{}}{\Sigma}$ for quadratic optimization, one need correct only the basis $\textsc{col}{\hskip 0.5pt(H)}$. This amounts to finding the (asymptotic) roots of the function $\scrE_{p}(\hskip 1.0pt\cdot\hskip 1.0pt)$. If the convergence to a root is sufficiently rapid, one may then estimate $\langle u_{H},\Gamma u_{H}\rangle$ closely by $\hat{\gamma}^{2}$ to bring the discrepancy $\hat{D}_{p}$ to one per Theorem 1. We conclude this section by showing that for many applications the rate of convergence of $\scrE_{p}(H)$ is less important than Theorem 1 suggests.

To illustrate some important examples in practice, we consider the following canonical, constrained optimization problem.

Now, $\hat{Q}(\hskip 1.0pt\cdot\hskip 1.0pt)$ in $(\ref{hQx})$ is the Lagrangian for $(\ref{minvar})$ with $c_{0}=0$ and $c_{1}=\langle\zeta,\bm{\hat{}}{\Sigma}^{-1}\zeta\rangle^{-1}$, which decays as $1/p$ under Assumption 2. The minimizer $\hat{w}\in\bbR^{p}$ of $(\ref{minvar})$ corresponds to the weights of a minimum variance portfolio of financial assets with $\zeta=1_{p}$ implementing the “full-investment” constraint. Minimum variance and the more general mean-variance optimized portfolios are widely used in finance. Here, the $p$ entries of a column of $B$ in $(\ref{bSig})$ represent the exposures of $p$ assets to that risk factor, e.g., market risk (bull/bear market), industry risk (energy, automotive, etc.), climate risk (migration, drought, etc.), innovation risk (Chat GPT, etc). Similar formulations based on $(\ref{minvar})$ arise in signal-noise processing and climate science (see Appendix E).¹⁴¹⁴14In signal-noise processing, $(\ref{minvar})$ maximizes the signal-to-noise ratio of a beamformer with $\zeta$ referred to as the “steering vector”. The same is done for optimal fingerprinting in climate science with $\zeta$ called the “guess pattern”. We review this literature with emphasis on estimation of $\Sigma$ in Appendix E.

Continuing with the above example, the minimum $\hat{\sigma}^{2}$ of $(\ref{minvar})$ corresponds to the variance of the estimated portfolio $\hat{w}$, while the expected out-of-sample variance is,

We have $\hat{D}_{p}=2-\hat{\mu}_{p}^{2}\scrV_{p}^{2}$ (see Appendix A) and, under the conditions of Theorem 1,

because $|\zeta|^{2}=|1_{p}|^{2}=p$. Because $|\Lambda_{p}|^{2}/p$ converges in $(0,\infty)$ as $p\uparrow\infty$ under our Assumption 1(b), we achieve (in expectation) an asymptotically riskless portfolio provided the convergence $|\scrE_{p}(H)|\to 0$ and irrespective of its rate.

We analyze the asymptotics $\scrE_{p}(H)$ for the pca estimate $H$. Lemma 2 with $K=\calS_{p}$ in $(\ref{pcs})$ and $\scrH^{\top}\scrH=I$ imply that $z_{H}=HH^{\dagger}z=\scrH\scrH^{\top}z=z_{\scrH}$ and $\langle z,z_{\scrH}\rangle=|z_{\scrH}|^{2}$, which reduces $(\ref{optbias})$ to

The unknowns in $(\ref{optbias_gps})$ are $\scrB^{\top}z$ and $\scrB^{\top}\scrH$ and we provide theoretical evidence that they cannot be estimated from data in Section 4 without very strong assumptions. Here, we nevertheless obtain an estimate of the length $|\scrE_{p}(H)|$. The following addresses a division by zero in $(\ref{optbias_gps})$. Recall that $z=\frac{\zeta}{|\zeta|}$ and $\zeta_{B}=BB^{\dagger}\zeta$ per $(\ref{eH})$.

Our next assumption concerns the matrices $\calX$ and $\calE$ in $(\ref{data})$. These guarantee that almost all realizations of the data $Y$ have full rank for sufficiently large $p$, allowing us to treat $n_{+}$ in $(\ref{snr-bulk})$ as $n_{+}=n$ when $J=I$ and $n_{+}=n-1$ otherwise.

These conditions are discussed below. Our fundamental result on pca (in conjunction with Theorem 1) may now be stated. Its proof is deferred to Appendix B.

We remark that $\varphi=\frac{\Pi\scrH^{\top}z}{|z-z_{\scrH}|}\in\bbR^{q}$ is computable solely from the data $Y$ with almost every $|\varphi|$ bounded in $[0,\infty)$ eventually. Theorem 3 demonstrates that pca, and sample eigenvectors $\scrH$ in particular, lead to poor “plug-in” covariance estimators for quadratic optimization unless every column of $\scrH$ is eventually orthogonal to $\zeta$. So typically, the discrepancy $\hat{D}_{p}$ in $(\ref{QhxD})$ between the estimated and realized optima diverges to $-\infty$ as $p$ grows and at a linear rate. In the portfolio application of Section 2.2, this covariance results in strictly positive expected (out-of-sample) portfolio risk $\scrV_{p}$ per $(\ref{tru})$–$(\ref{tru_asymp})$ asymptotically, which may be approximated by using $|\varphi|$.

We make some remarks on Assumption 6. Conditions (a)–(c) are straightforward, but we mention that the invertibility of $\calX^{\top}J\calX$ is closely related to the requirement that $n_{+}>q$ in condition (b). Condition (c) fails when $g$ in $(\ref{Jc})$ lies in $\textsc{col}{\hskip 0.5pt(\calX)}$ but such a case is dealt with by rewriting the data in $(\ref{data})$ as $Y=\alpha g^{\top}+B_{\alpha}\calX^{\top}_{\alpha}+\calE$ for some $B_{\alpha}$ and $\calX_{\alpha}$ of $q-1$ columns each, and some mean vector $\alpha\in\bbR^{p}$. Then, we have $YJ=B_{\alpha}\calX^{\top}_{\alpha}J+\calE J$, and it only remains to check if condition (c) holds with the matrix $\calX_{\alpha}$ replacing $\calX$. Conditions (d)–(f) require that strong laws of large numbers hold for the columns of the sequence $\calE=\calE_{p\times n}$. These roughly state that the columns of $\calE$ are stationary with weakly dependent entries having bounded fourth moments. All three are easily verified for the $\calE_{p\times n}$ populated by i.i.d. Gaussian random entries. Lastly, we remark that if conditions (e) and (f) hold for $J=I$ they hold for any $J$.