MNRASDate MNRASYear MNRASDate

# A comparison of shrinkage estimators of the cosmological precision matrix

Marnix J. Looijmans, 1Mike (Shengbo) Wang 1 and Florian Beutler 1
1Institute for Astronomy, University of Edinburgh, Royal Observatory
Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK Email:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form 21 February 2024)

###### Abstract

The determination of the covariance matrix and its inverse, the _precision_
matrix, is critical in the statistical analysis of cosmological measurements.
The covariance matrix is typically estimated with a limited number of
simulations at great computational cost before inversion into the precision
matrix; therefore, it can be ill-conditioned and overly noisy when the sample
size $n$ used for estimation is not much larger than the data vector
dimension. In this work, we consider a class of methods known as _shrinkage_
estimation for the precision matrix, which combines an empirical estimate with
a target that is either analytical or stochastic. These methods include linear
and non-linear shrinkage applied to the covariance matrix (the latter
represented by the so-called NERCOME estimator), and the direct linear
shrinkage estimation of the precision matrix which we introduce in a
cosmological setting. Using Bayesian parameter inference as well as metrics
like matrix loss functions and the eigenvalue spectrum, we compare their
performance against the standard sample estimator with varying sample size
$n$. We have found the shrinkage estimators to significantly improve the
posterior distribution at low $n$, especially for the linear shrinkage
estimators either inverted from the covariance matrix or applied directly to
the precision matrix, with an empirical target constructed from the sample
estimate. Our results should be particularly relevant to the analyses of
Stage-IV spectroscopic galaxy surveys such as the Dark Energy Spectroscopic
Instrument (DESI) and Euclid, whose statistical power can be limited by the
computational cost of obtaining an accurate precision matrix estimate.

###### keywords:

methods: statistical – methods: data analysis – methods: numerical – large-
scale structure of the Universe

††pubyear: MNRASYear21 February 2024††pagerange: A comparison of shrinkage
estimators of the cosmological precision matrix–REFERENCES

##  1 Introduction

With the arrival of Stage-IV spectroscopic galaxy surveys such as the Dark
Energy Spectroscopic Instrument (DESI)111desi.lbl.gov (DESI Collaboration,
2024, 2023) and Euclid222sci.esa.int/euclid, euclid-ec.org (Laureijs et al.,
2011), we will be able to observe cosmic structure on scales larger than ever
before. At the same time, there has been significant development in
perturbative models of non-linear clustering on smaller scales (e.g. Perko et
al., 2016; d’Amico et al., 2020; Semenaite et al., 2022). The wider dynamic
range accessible to galaxy surveys poses a statistical challenge as the number
of modes of cosmological fluctuations becomes very large, especially for the
determination of the associated covariance matrix that is vital to most data
analyses.

The covariance matrix characterizes the uncertainties and thus the
constraining power of the analysis. Whilst it is sometimes possible to derive
an analytical expression from perturbation theory (e.g. Grieb et al., 2016; Li
et al., 2019; Sugiyama et al., 2020; Wadekar et al., 2020), including higher-
order corrections and accounting for systematic effects such as the survey
window is far from straightforward (Wadekar & Scoccimarro, 2020). The
alternative is an empirical approach whereby a sample estimate of the
covariance matrix is obtained from a limited number of high-fidelity mock
realizations. The main challenge with this method is that the number of
realizations depends on the dimension of the data vector, which can mean that
thousands of mock catalogues are needed. Because the evaluation of the
likelihood function typically requires the inversion of the covariance matrix
into the precision matrix, the covariance matrix estimate must not be
singular, which is the case when the number of realizations is no greater than
the data vector dimension. Previous studies (Taylor et al., 2013; Taylor &
Joachimi, 2014; Sellentin & Heavens, 2016) have also shown that using only a
limited number of data realizations to estimate the covariance matrix induces
noise to which the operation of matrix inversion is highly sensitive.
Moreover, since matrix inversion is a non-linear operation, an unbiased
estimator of the covariance matrix does not generally result in an unbiased
estimator of the precision matrix (Hartlap et al., 2006).

To alleviate the issues associated with covariance matrix estimation,
resampling methods such as jackknife and bootstrapping have been proposed
(e.g. Friedrich et al., 2015; Escoffier et al., 2016; O’Connell & Eisenstein,
2019; Philcox et al., 2019; Mohammad & Percival, 2022), which construct the
covariance matrix from sub-samples of the observed data itself. However, these
sub-samples are obtained by dividing the original simulations into much
smaller sub-volumes, which do not necessarily contain the largest modes of
interest.

In this work, we investigate an alternative class of methods known as
_shrinkage_ estimation, which combines empirical estimates of the covariance
matrix with a _target_ that is often (but not always) an analytical proxy. The
standard linear shrinkage method in the statistical literature (Ledoit & Wolf,
2003, 2004; Schäfer & Strimmer, 2005) has been previously explored in a
cosmological context by Pope & Szapudi (2008), and the NERCOME333Non-
parametric Eigenvalue-Regularized COvariance Matrix Estimator estimator from
non-linear shrinkage methods (Ledoit & Wolf, 2012) has been first introduced
by Joachimi (2017) in the same context. In both methods, the shrinkage
estimate of the covariance matrix still needs to be inverted into the
precision matrix estimate for likelihood analysis, and the effect of inversion
is not yet fully understood. We will compare the performance of these two
estimation methods with a third method, introduced and applied here in a
cosmological context for the first time, which applies linear shrinkage
directly to precision matrix estimation instead (Bodnar et al., 2016).

This paper is organized as follows: we introduce the different covariance and
precision matrix estimators in section 2, and apply them to a power spectrum
analysis in section 3; we then compare their performance in section 4, and
discuss our findings before concluding in section 5.

For reference and clarity, our notation is listed in Table 1.

Table 1: Notation used in this work. Symbol | Definition
---|---
$\operatorname{diag}{\mathsfbf{A}}$ | Diagonal matrix consisting of the diagonal elements of $\mathsfbf{A}$
$\|\mathsfbf{A}\|_{\mathrm{F}}$ | Frobenius norm given by $\mathopen{}\mathclose{{}\left[\operatorname{tr}\\!{\mathopen{}\mathclose{{}\left(\mathsfbf{A}{\mathsfbf{A}}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\right)}}\right]^{1/2}$
$\|\mathsfbf{A}\|_{\mathrm{tr}}$ | Trace norm given by $\operatorname{tr}\\!{\bigl{[}\mathopen{}\mathclose{{}\left(\mathsfbf{A}{\mathsfbf{A}}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\right)^{1/2}\bigr{]}}$
$\mathopen{}\mathclose{{}\left\llbracket a,b}\right\rrbracket$ | Closed integer interval
$d$ | Length of the data vector
$n$ | Number of data realizations
$p$ | Number of inferred parameters
$\mathsfbf{X}\in{}^{d\times n}$ | Data matrix
$\mathsfbf{\Sigma}\in{}^{d\times d}$ | True covariance matrix
$\mathsfbf{\Psi}\in{}^{d\times d}$ | True precision matrix given by $\mathsfbf{\Sigma}^{-1}$
$\mathsfbf{S}$ | Sample estimator of $\mathsfbf{\Sigma}$
$\mathsfbf{\Pi}$ | Sample estimator of $\mathsfbf{\Psi}$
$\mathsfbf{T}$ | Covariance matrix target
$\mathsfbf{\Pi}_{\mkern 2.0mu0}$ | Precision matrix target
$\hat{\mathsfbf{\Sigma}}_{\mathrm{LS}}$ | Linear shrinkage estimator of $\mathsfbf{\Sigma}$
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ |  $\hat{\mathsfbf{\Sigma}}_{\mathrm{LS}}^{-1}$, inverted covariance shrinkage estimator of $\mathsfbf{\Psi}$
$\hat{\mathsfbf{\Sigma}}_{\mathrm{NLS}}$ | Non-linear shrinkage (NERCOME) estimator of $\mathsfbf{\Sigma}$
$\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$ |  $\hat{\mathsfbf{\Sigma}}_{\mathrm{NLS}}^{-1}$, NERCOME estimator of $\mathsfbf{\Psi}$
$\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}$ | Direct linear shrinkage estimator of $\mathsfbf{\Psi}$
$\lambda$ | Linear-shrinkage intensity parameter
$\alpha,\beta$ | Linear-shrinkage mixing parameters
$s$ | NERCOME split parameter

For reproducibility and reference, the implementation of the different
estimation methods considered in this work can be found in our public
repository online (see Data Availability for details).

##  2 Estimation Methods

In this section, we describe the different precision matrix estimation methods
to be compared. Most of them invert a covariance matrix estimate at the final
stage, but the new method we adopt involves matrix inversion only as an
intermediate step. The implications of covariance matrix inversion will be
discussed.

###  2.1 Sample estimation

Suppose we have a data matrix $\mathsfbf{X}\in{}^{d\times n}$ consisting of
$n$ independent realizations of a random vector of length $d$, which has an
underlying covariance matrix $\mathsfbf{\Sigma}\in{}^{d\times d}$ that is a
priori unknown. Without loss of generality, we assume the data vectors are
mean-subtracted. The standard sample estimator of the covariance matrix is
then given by

$\mathsfbf{S}=\frac{1}{n-1}\mathsfbf{X}{\mathsfbf{X}}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\,.$
(1)

It is unbiased with expectation
$\operatorname{\mathbb{E}}\\!{\mathopen{}\mathclose{{}\left[\mathsfbf{S}}\right]}=\mathsfbf{\Sigma}$.

In a typical likelihood function, it is the precision matrix that appears
instead of the covariance matrix, so one may wish to estimate the true
precision matrix $\mathsfbf{\Psi}\equiv\mathsfbf{\Sigma}^{-1}$. From equation
(1), we obtain the sample precision matrix estimator

$\mathsfbf{\Pi}=\gamma\mathsfbf{S}^{-1}\,,$ (2)

where the prefactor

$\gamma=\frac{n-d-2}{n-1}$ (3)

corrects for the multiplicative bias resulting from the inversion of the
covariance matrix estimator $\mathsfbf{S}$, which is a non-linear operation.
This prefactor $\gamma$ is commonly known as the ‘Hartlap factor’ in
cosmological analyses (Hartlap et al., 2006). It is immediately clear that the
number of data realizations must satisfy $n>d+2$, and in the limit
$n\to\infty$, $\mathsfbf{S}^{-1}$ represents an asymptotically unbiased
precision matrix estimator in itself.

###  2.2 Linear shrinkage

The _linear shrinkage_ (LS) method uses a convex combination of an empirical
covariance matrix estimate $\hat{\mathsfbf{\Sigma}}$ and a _target_ matrix
$\mathsfbf{T}$ to construct a new covariance matrix estimator (Ledoit & Wolf,
2003, 2004),

$\hat{\mathsfbf{\Sigma}}_{\mathrm{LS}}=(1-\lambda)\hat{\mathsfbf{\Sigma}}+\lambda\mathsfbf{T}\,,$
(4)

where $\lambda\in[0,1]$ is the shrinkage _intensity_. The empirical estimate
$\hat{\mathsfbf{\Sigma}}$ is commonly chosen to be the _unbiased_ sample
covariance matrix estimate $\mathsfbf{S}$ in equation (1), which we shall
assume for the rest of this work; the target $\mathsfbf{T}$ is often, but not
always, an analytical quantity.

The idea behind shrinkage methods is that the target matrix has a smaller
variance (possibly zero) than the empirical estimate, but may be otherwise
biased, and one aims to find the optimal trade-off between bias and variance.
Following Pope & Szapudi (2008), we first define

$\widebar{W}_{ij}=\frac{1}{n}\sum_{k=1}^{n}W_{ij}^{(k)}\quad\text{with}\quad
W_{ij}^{(k)}=x_{i}^{(k)}x_{j}^{(k)}\,,$ (5)

where $x_{i}^{(k)}$ is the $(i,k)$-th element of $\mathsfbf{X}$. The elements
of the sample covariance matrix are then given by

$S_{ij}=\frac{n}{n-1}\widebar{W}_{ij}\,,$ (6)

whose covariance can be estimated by

$\operatorname{\widehat{\mathrm{Cov}}}\\!{\mathopen{}\mathclose{{}\left(S_{ij},S_{lm}}\right)}=\frac{n}{(n-1)^{3}}\sum_{k=1}^{n}\mathopen{}\mathclose{{}\left(W_{ij}^{(k)}-\widebar{W}_{ij}}\right)\mathopen{}\mathclose{{}\left(W_{lm}^{(k)}-\widebar{W}_{lm}}\right)\,.$
(7)

Ledoit & Wolf (2003) and Schäfer & Strimmer (2005) have shown that the optimal
intensity parameter is estimated by

$\hat{\lambda}^{*}=\frac{\sum_{i,j}\mathopen{}\mathclose{{}\left[\operatorname{\widehat{\mathrm{Var}}}\\!{\mathopen{}\mathclose{{}\left(S_{ij}}\right)}-\operatorname{\widehat{\mathrm{Cov}}}\\!{\mathopen{}\mathclose{{}\left(S_{ij},T_{ij}}\right)}}\right]}{\sum_{i,j}\mathopen{}\mathclose{{}\left({S}_{ij}-T_{ij}}\right)^{2}}\,,$
(8)

where
$\operatorname{\widehat{\mathrm{Var}}}\\!{\mathopen{}\mathclose{{}\left(S_{ij}}\right)}\equiv\operatorname{\widehat{\mathrm{Cov}}}\\!{\mathopen{}\mathclose{{}\left(S_{ij},S_{ij}}\right)}$.
If $\hat{\lambda}^{*}$ is found to be negative, it is then clipped to
$\hat{\lambda}^{*}=0$; if it is greater than 1, it is clipped to
$\hat{\lambda}^{*}=1$.

Finally, an estimate of the precision matrix is obtained by inversion,
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}=\hat{\mathsfbf{\Sigma}}_{\mathrm{LS}}^{-1}$.
Since there is no general analytical form for the probability distribution of
$\hat{\mathsfbf{\Sigma}}_{\mathrm{LS}}$, unlike for the sample estimator
$\mathsfbf{S}$ which follows the Wishart distribution (Hartlap et al., 2006),
we do not include an analogous Hartlap factor here (de la Torre et al., 2013).

###  2.3 Non-linear shrinkage

The linear shrinkage estimator is in fact a first-order approximation to a
more general class of _non-linear shrinkage_ (NLS) estimators that attempt to
shrink the dynamic range of the eigenvalues of the covariance matrix estimate
and thus mitigate ill-conditioning when the data vector dimension becomes
large (Ledoit & Wolf, 2012). In this work, we consider in particular the
_NERCOME estimator_ of the covariance matrix proposed by Ledoit & Wolf (2012),
Abadir et al. (2014) and Lam (2016), which has been previously applied by
Joachimi (2017) and Gouyou Beauchamps et al. (2023) in a cosmological setting.

The NERCOME estimator is constructed from a set of data realizations as
follows:

  1. 1.

Without necessarily preserving their order, the columns of the data matrix
$\mathsfbf{X}$ are split into a $d\times s$ matrix $\mathsfbf{X}_{1}$ and a
$d\times(n-s)$ matrix $\mathsfbf{X}_{2}$ for a given value of the split
parameter $s\in\mathopen{}\mathclose{{}\left\llbracket
2,n-2}\right\rrbracket$;

  2. 2.

Two sample estimates of the covariance matrix, denoted $\mathsfbf{S}_{1}$ and
$\mathsfbf{S}_{2}$, are obtained from $\mathsfbf{X}_{1}$ and
$\mathsfbf{X}_{2}$ respectively;

  3. 3.

The sample estimates are diagonalized,
$\mathsfbf{S}_{a}=\mathsfbf{U}_{a}\mathsfbf{D}_{a}{\mathsfbf{U}}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}_{a}$
with $a\in\\{1,2\\}$, where $\mathsfbf{D}_{a}$ is the diagonal matrix of
eigenvalues and $\mathsfbf{U}_{a}$ is the orthogonal matrix of eigenvectors;

  4. 4.

The quantity
$\mathsfbf{Z}=\mathsfbf{U}_{1}\operatorname{diag}\\!{\mathopen{}\mathclose{{}\left({\mathsfbf{U}}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}_{1}\mathsfbf{S}_{2}\mathsfbf{U}_{1}}\right)}{\mathsfbf{U}}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}_{1}$
is computed for a given split at the chosen split parameter value $s$;

  5. 5.

Over all $\binom{n}{s}$ possible choices of the split into
$\mathsfbf{X}_{1,2}$ for the chosen value of $s$, the average of
$\mathsfbf{Z}$ is calculated to be $\widebar{\mathsfbf{Z}}(s)$, and
correspondingly the sample estimates $\mathsfbf{S}_{2}$ are averaged to obtain
$\widebar{\mathsfbf{S}}_{2}(s)$;

  6. 6.

The _distance_ function

$Q(s)=\big{\|}\widebar{\mathsfbf{Z}}(s)-\widebar{\mathsfbf{S}}_{2}(s)\big{\|}_{\mathrm{F}}$
(9)

(or equivalently $Q^{2}$) is minimized to find an optimal split value $s^{*}$.

The NERCOME estimator of the covariance matrix is then defined as
$\hat{\mathsfbf{\Sigma}}_{\mathrm{NLS}}=\widebar{\mathsfbf{Z}}(s^{*})$. As
before, the precision matrix estimator is given by its inverse,
$\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}=\widebar{\mathsfbf{Z}}(s^{*})^{-1}$, and
we do not include a Hartlap-like factor.

###  2.4 Direct linear shrinkage for the precision matrix

Instead of inverting the linear shrinkage estimator of the covariance matrix,
one can apply shrinkage directly to an empirical estimate and a target for the
precision matrix. Here we consider the inverse of the sample covariance matrix
estimate, $\mathsfbf{S}^{-1}$, as the empirical estimate, with some suitable
precision matrix target denoted by $\mathsfbf{\Pi}_{\mkern 2.0mu0}$ to give
the direct linear shrinkage estimator of the precision matrix:

$\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}=\alpha\mathsfbf{S}^{-1}+\beta\mkern
1.0mu\mathsfbf{\Pi}_{\mkern 2.0mu0}\,.$ (10)

Here $\alpha$ and $\beta$ are the _mixing parameters_ , and $n>d$ must be
satisfied to ensure that $\mathsfbf{S}$ is non-singular; in addition, the
target $\mathsfbf{\Pi}_{\mkern 2.0mu0}$ is required to satisfy the condition
$\sup_{d}d^{-1}\|\mathsfbf{\Pi}_{\mkern 2.0mu0}\|_{\mathrm{F}}^{2}\leqslant a$
for some constant $a>0$. Bodnar et al. (2016) have shown that the optimal
mixing parameters can be estimated by

$\displaystyle\hat{\alpha}^{*}$
$\displaystyle=1-\frac{d}{n}-\frac{n^{-1}\|\mathsfbf{S}^{-1}\|_{\mathrm{tr}}^{2}\mkern
2.0mu\|\mathsfbf{\Pi}_{\mkern
2.0mu0}\|_{\mathrm{F}}^{2}}{\|\mathsfbf{S}^{-1}\|_{\mathrm{tr}}^{2}\mkern
2.0mu\|\mathsfbf{\Pi}_{\mkern
2.0mu0}\|_{\mathrm{F}}^{2}-\mathopen{}\mathclose{{}\left[\operatorname{tr}\\!{\mathopen{}\mathclose{{}\left(\mathsfbf{S}^{-1}\mathsfbf{\Pi}_{\mkern
2.0mu0}}\right)}}\right]^{2}}\,,$ (11a) $\displaystyle\hat{\beta}^{*}$
$\displaystyle=\mathopen{}\mathclose{{}\left(1-\frac{d}{n}-\hat{\alpha}^{*}}\right)\frac{\operatorname{tr}\\!{\bigl{(}\mathsfbf{S}^{-1}\mathsfbf{\Pi}_{\mkern
2.0mu0}\bigr{)}}missing}{\|\mathsfbf{\Pi}_{\mkern
2.0mu0}\|_{\mathrm{F}}^{2}}\,,$ (11b)

and the estimator $\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}$ improves
asymptotically as $n\to\infty$ with $\ifrac{d}{n}\to c\in(0,1)$ for some fixed
constant $c$.

##  3 Application to Power Spectrum Analysis

To demonstrate the performance of the estimation methods reviewed in the
previous section, we apply them to a galaxy clustering power spectrum analysis
of the Baryon Oscillations Spectroscopic Survey (BOSS) Data Release 12 (DR12)
catalogues (BOSS Collaboration et al., 2015; Reid et al., 2016) for the North
Galactic Cap (NGC) in the redshift bin at $z=0.38$, with the accompanying
$2048$ Patchy mock catalogues (Kitaura et al., 2016) for covariance matrix
estimation. The data set we use includes measurements of the clustering power
spectrum monopole $P_{0}$ and quadrupole $P_{2}$ over the wavenumber range
$k\in[0.01,0.1]$\mathrm{\text{$h$}}\text{\,}{\mathrm{Mpc}}^{-1}$$ in bins of
width
$\Delta{k}=$0.01\text{\,}\mathrm{\text{$h$}}\text{\,}{\mathrm{Mpc}}^{-1}$$,
which results in data vectors of length $d=18$.444Full data set including the
transformation matrices are available at fbeutler.github.io/hub/deconv_paper.

Given the wavenumber range considered in this work, we adopt a linear
redshift-space power spectrum model (Kaiser, 1987) computed from Eisenstein &
Hu (1998), with the fiducial cosmology taken from the Planck 2015 results
(Planck Collaboration, 2016). The power spectrum model includes the monopole,
quadrupole and hexadecapole given by

$\displaystyle P_{0}(k)$
$\displaystyle=\mathopen{}\mathclose{{}\left(b^{2}+\frac{2}{3}bf+\frac{1}{5}f^{2}}\right)P_{\mathrm{m}}(k)\,,$
(12a) $\displaystyle P_{2}(k)$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{4}{3}bf+\frac{4}{7}f^{2}}\right)P_{\mathrm{m}}(k)\,,$
(12b) $\displaystyle P_{4}(k)$
$\displaystyle=\frac{8}{35}f^{2}P_{\mathrm{m}}(k)\,,$ (12c)

where $P_{\mathrm{m}}$ is the linear matter power spectrum, $b$ is the linear
bias and $f$ is the linear growth rate. We include the wide-angle and survey
window effects as matrices $\mathsfbf{W}$ and $\mathsfbf{M}$ respectively
(Beutler & McDonald, 2021), which transform the power spectrum model vector
$\textbf{{P}}_{\mathrm{model}}\mapsto\mathsfbf{W}\mathsfbf{M}\textbf{{P}}_{\mathrm{model}}$
evaluated in the same bins as for the measurements.

To assess the performance of the different estimators with varying sample
size, we have created a superset of $2040$ power spectrum measurements from
the $2048$ mock catalogues available. This superset is divided into $85$
subsets of $24$ mock measurements (i.e. $85$ data matrices
$\mathsfbf{X}\in{}^{18\times 24}$), or $68$ subsets of $30$ mock measurements
each (i.e. $68$ data matrices $\mathsfbf{X}\in{}^{18\times 30}$), so that we
can test covariance and precision matrix estimation with $n=$24$$ and $n=$30$$
sample sizes respectively.

###  3.1 Sample estimates

For each subset of the mock measurements, we compute a sample covariance
matrix estimate $\mathsfbf{S}$ using equation (1), which is then inverted into
a sample precision matrix estimate $\mathsfbf{\Pi}$ with the Hartlap factor
$\gamma$ being $0.174$ for $n=$24$$ and $0.344$ for $n=$30$$ respectively.

As a _reference_ covariance matrix, we also compute the sample covariance
matrix estimate $\mathsfbf{S}_{\mathrm{ref}}$ from all $2048$ mock
measurements, which is used to benchmark all other covariance matrix
estimators. The corresponding sample precision matrix estimate
$\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{ref}}$ is used as a reference precision
matrix. This has a Hartlap factor $\gamma=$0.991$\simeq 1$ as $n=$2048$\gg
d=18$, and we regard it as a good proxy for the unknown true precision matrix
$\mathsfbf{\Psi}$.

###  3.2 Inverted covariance shrinkage estimates

For the linear shrinkage estimator (4) of the covariance matrix, we consider
two possibilities for the target matrix $\mathsfbf{T}$. The first choice,
proposed by Schäfer & Strimmer (2005), is simply the diagonal of the sample
covariance matrix estimate:

$\mathsfbf{T}^{(1)}=\operatorname{diag}{\mathsfbf{S}}\,.$ (13)

In this case, only summands with $i\neq j$ contribute to the numerator and
denominator of expression (8). The second choice of target, following Hamilton
et al. (2006) and Pope & Szapudi (2008), is given by

$\big{(}\mathsfbf{T}^{(2)}\big{)}_{ij}=\frac{2}{N_{i}}[P_{\mathrm{fid}}(k_{i})]^{2}\delta^{\mathrm{(K)}}_{ij}\,,$
(14)

where $P_{\mathrm{fid}}$ is a fiducial power spectrum model for which we set
$b=2$ and $f=0.7$ and transform using the wide-angle and survey window
matrices, and $N_{i}$ is the number of modes in the bin at effective
wavenumber $k_{i}$. Since $\mathsfbf{T}^{(2)}$ is deterministic,
$\operatorname{\widehat{\mathrm{Cov}}}\\!{\mathopen{}\mathclose{{}\left(S_{ij},T_{ij}}\right)}$
identically vanishes in equation (8).

For both choices of the target matrix, we compute the linear shrinkage
estimate $\hat{\mathsfbf{\Sigma}}_{\mathrm{LS}}$ of the covariance matrix from
each subset of mock measurements; this gives $85$ estimates with $n=$24$$ and
$68$ estimates with $n=$30$$. The distributions of the optimal shrinkage
intensity parameter $\hat{\lambda}^{*}$ for these two target choices are
compared in Fig. 1.

Figure 1:  Distributions of the optimal shrinkage intensity
$\hat{\lambda}^{*}$ for target matrix choices $\mathsfbf{T}^{(1)}$ versus
$\mathsfbf{T}^{(2)}$. The top panel shows the $85$ intensity parameters
computed from subsets of $n=24$ mock measurements; the bottom panel shows the
$68$ intensity parameters computed for $n=30$. The diagonal dashed lines show
the line of equality for $\hat{\lambda}^{*}$ between the two target choices.

We have found that the shrinkage intensity $\hat{\lambda}^{*}$ with target
$\mathsfbf{T}^{(1)}$ is equal to unity for $31$ out of $85$ data matrices with
$n=$24$$ and for $21$ out of $68$ data matrices with $n=$30$$. This indicates
that in these cases, the linear shrinkage estimate in equation (4) is
dominated by the target with no contribution from the sample estimate. In
contrast, with target $\mathsfbf{T}^{(2)}$, the shrinkage intensity
$\hat{\lambda}^{*}$ favours smaller values, which means that the target tends
to contribute less to the shrinkage estimate than the sample estimate does. If
the reference covariance matrix $\mathsfbf{S}_{\mathrm{ref}}$ is used as
$\mathsfbf{S}$ in the linear shrinkage estimator, the optimal intensity
parameter is $\hat{\lambda}^{*}\approx$0.028$$ for target $\mathsfbf{T}^{(1)}$
and $\hat{\lambda}^{*}\approx$0.002$$ for target $\mathsfbf{T}^{(2)}$. In both
cases, the linear shrinkage estimator is almost fully determined by the sample
estimate of the covariance matrix. This is expected since our reference
covariance matrix $\mathsfbf{S}_{\mathrm{ref}}$ is supposed to closely
approximate the true covariance matrix $\mathsfbf{\Sigma}$.

To obtain the corresponding precision matrix estimates
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$, the aforementioned linear-shrinkage
covariance matrix estimates are simply inverted.

###  3.3 NERCOME estimates

We follow the steps detailed in section 2.3 to compute $85$ NERCOME covariance
matrix estimates $\hat{\mathsfbf{\Sigma}}_{\mathrm{NLS}}$ from the data
matrices $\mathsfbf{X}\in{}^{18\times 24}$ and $68$ estimates from the data
matrices $\mathsfbf{X}\in{}^{18\times 30}$. For reference, we also compute a
single NERCOME estimate using all $2048$ mock measurements. As before, the
corresponding precision matrix estimates
$\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$ are obtained by matrix inversion.

A caveat here is that since $\binom{n}{s}$ can be very large for some values
of $s$, we only randomly draw
$N_{\mathrm{draw}}=\min\mathopen{}\mathclose{{}\left\\{\binom{n}{s},$1000$}\right\\}$
splits of the data matrix $\mathsfbf{X}$ at step (v) in section 2.3 to compute
$\widebar{\mathsfbf{Z}}(s)$ and $\widebar{\mathsfbf{S}}_{2}(s)$ (Joachimi,
2017).

To minimize the distance function $Q(s)$ in equation (9), we evaluate it at
each point in the integer interval
$s\in\mathopen{}\mathclose{{}\left\llbracket 2,n-2}\right\rrbracket$, except
when $n=2048$, in which case we evaluate $Q(s)$ only at
$s\in\mathopen{}\mathclose{{}\left\\{0.1n,0.15n,0.2n,\dots,0.9n,2\sqrt{n},n-1.5\sqrt{n},\,n-2.5\sqrt{n}}\right\\}$,
as proposed by Joachimi (2017) and Lam (2016).

The distribution of the optimal split parameters $s^{*}$ for all NERCOME
estimates computed with $n=$24$$ and $n=$30$$ is shown in Fig. 2.

Figure 2:  Normalized histogram showing the distribution of the optimal split
parameter $s^{*}$ that minimizes the distance function $Q(s)$ (equation 9) as
a ratio $\ifrac{s}{n}$ to the number of data realizations. The top panel
corresponds to $85$ NERCOME covariance matrix estimates computed from subsets
of $n=$24$$ mock measurements; the bottom panel corresponds to $68$ estimates
computed for $n=$30$$. The vertical dashed line shows the optimal split
parameter for $n=$2048$$ mock measurements.

The split-parameter distribution indicates that $\ifrac{s^{*}}{n}\simeq 0.8$
is favoured, which is slightly different to the fixed split
$\ifrac{s^{*}}{n}=\ifrac{2}{3}$ found in Joachimi (2017). For reference, the
split parameter for the NERCOME estimate computed using all $2048$ mocks is
$\ifrac{s^{*}}{n}=\ifrac{$1638$}{$2048$}\simeq$0.80$$.

###  3.4 Direct precision shrinkage estimates

Finally, for the direct linear shrinkage estimator (10) of the precision
matrix, we consider two possibilities for the target matrix
$\mathsfbf{\Pi}_{\mkern 2.0mu0}$. The first choice is the matrix consisting of
eigenvalues of the inverse sample covariance matrix estimate on the diagonal,
i.e.

$\big{(}\mathsfbf{\Pi}_{\mkern
2.0mu0}^{(1)}\big{)}_{ij}=e_{i}\delta^{\mathrm{(K)}}_{ij},$ (15)

where $\\{e_{i}\\}_{i=1}^{d}$ are the eigenvalues of $\mathsfbf{S}^{-1}$ in
ascending order. The second choice is to invert the analytical covariance
matrix target from equation (14),

$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}=\big{(}\mathsfbf{T}^{(2)}\big{)}^{-1}.$
(16)

We then compute two linear shrinkage estimates $\mathsfbf{\Pi}_{\mkern
2.0mu\mathrm{LS}}$ of the precision matrix, one per target, for each subset of
mock measurements with $n=$24$$ or $n=$30$$ using equation (10).

The estimated optimal shrinkage mixing parameters $\hat{\alpha}^{*}$ and
$\hat{\beta}^{*}$ are computed using equations (11) and shown in Fig. 3.

Figure 3:  Distributions of the optimal linear-shrinkage mixing parameters
$\hat{\alpha}^{*}$ and $\hat{\beta}^{*}$ for two choices of the target matrix
$\mathsfbf{\Pi}_{\mkern 2.0mu0}$ used in the direct precision matrix estimate
$\hat{\mathsfbf{\Pi}}_{\mathrm{LS}}$. Each top panel shows the $85$ mixing
parameters computed for subsets of $n=$24$$ mock measurements, and each bottom
panel shows the $68$ mixing parameters computed for $n=$30$$. Note the
different scales on the vertical axes.

We note that $\hat{\alpha}^{*}$ and $\hat{\beta}^{*}$ have an approximate
linear relationship when $\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(1)}$ is used as
target. When we use $\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}$ as target, the
linear relationship disappears. For target $\mathsfbf{\Pi}_{\mkern
2.0mu0}^{(1)}$, we have noted that $\hat{\alpha}^{*}$ is negative for one out
of the $85$ direct precision shrinkage estimates with $n=$24$$, indicating
that the signs of all the entries in $\mathsfbf{S}^{-1}$ are flipped for this
estimate. The mixing parameter $\hat{\alpha}^{*}$ is never negative for any
other number of mock measurements $n$ or target $\mathsfbf{\Pi}_{\mkern
2.0mu0}^{(2)}$. Moreover, it can be seen from the distribution of mixing
parameters that $\hat{\alpha}^{*}$ becomes larger when a larger number of mock
measurements is used. This indicates that the inverse sample covariance matrix
estimate $\mathsfbf{S}^{-1}$ has greater weight in shrinkage with larger
numbers of mock measurements, as one would expect. Finally, the values of
$\hat{\beta}^{*}$ for the target $\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}$ are
smaller by an order of magnitude than those for the the target
$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(1)}$, implying that the analytical target is
generally less favoured than the empirical target.

For reference, we have also listed in Table 2 the optimal shrinkage mixing
parameters when $\mathsfbf{S}_{\mathrm{ref}}$ is used for the precision matrix
estimate in equation (10). As expected, the $\hat{\alpha}^{*}$ values are
close to 1 and the $\hat{\beta}^{*}$ values are close to 0.

Table 2:  Linear-shrinkage mixing parameters for the direct precision matrix estimate computed with $\mathsfbf{S}=\mathsfbf{S}_{\mathrm{ref}}$ in equation (10) for different target matrices.  Target | $\hat{\alpha}^{*}$ |  $\hat{\beta}^{*}$
---|---|---
$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(1)}$ (empirical) | $0.989$ |  $4.71\text{\times}{10}^{-4}$
$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}$ (analytical) | $0.989$ |  $1.39\text{\times}{10}^{-5}$

##  4 Performance Comparison

We will next compare the performance of seven types of precision matrix
estimate constructed from the previous section for application to the BOSS
DR12 power spectrum analysis, which are listed below for clarity and
completeness:

  1. 1.

the sample estimate $\mathsfbf{\Pi}$ introduced in section 2.1 and computed in
section 3.1;

  2. 2.

the inverted covariance linear-shrinkage estimate,
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$, introduced in section 2.2 with the
empirical target $\mathsfbf{T}^{(1)}$ or the analytical target
$\mathsfbf{T}^{(2)}$ defined in section 3.2 (hereafter referred to as simply
the covariance shrinkage estimate);

  3. 3.

the inverted NERCOME estimate, $\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$,
introduced in section 2.3 and constructed in section 3.3;

  4. 4.

the direct precision linear-shrinkage estimate, $\mathsfbf{\Pi}_{\mkern
2.0mu\mathrm{LS}}$, introduced in section 2.4 with the empirical target
$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(1)}$ or the analytical target
$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}$ defined in section 3.4 (hereafter
referred to as simply the precision shrinkage estimate).

To this end, two types of performance metrics are considered: the first is
based on the eigenvalue spectrum and loss function of the matrices themselves;
the second is based on parameter inference, which is the ultimate objective of
precision matrix estimation in cosmological likelihood analysis. The benchmark
used in all comparisons is the ‘reference’ precision matrix
$\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{ref}}$, which we recall from section 3 is
the sample estimate computed using all $2048$ mock measurements and treated as
a proxy for the true precision matrix $\mathsfbf{\Psi}$.

###  4.1 Eigenvalue spectrum and loss function

In general, since covariance and precision matrices are inversely related,
smaller eigenvalues of the precision matrix correspond to larger inferred
uncertainties, and the lowest eigenvalues correspond to the highest-variance
principal components (eigenvectors) of the measurements. The eigenvalue
spectra of the different precision matrix estimates computed from subsets of
$n=$24$$ mock measurements are shown in Fig. 4, and those from subsets of
$n=$30$$ mock measurements are shown in Fig. 5. For $n=$2048$$ where all mock
measurements are used in constructing the estimates (i.e. computed from the
reference covariance matrix $\mathsfbf{S}_{\mathrm{ref}}$ wherever
applicable), we show in Fig. 6 the relative differences between the eigenvalue
spectra instead.

Figure 4:  Eigenvalue spectra of the precision matrix estimates computed from
subsets of $n=$24$$ mock measurements. Each sub-figure shows the scatter of
$85$ spectra of one type of precision matrix estimate (with the target choice
indicated where applicable), which are compared with the spectrum of the
reference precision matrix $\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{ref}}$ shown
by the dashed red line.  Figure 5:  The same as Fig. 4, but with $68$
estimates computed from subsets of $n=$30$$ mock measurements in each sub-
figure.  Figure 6:  Relative difference,
$\ifrac{e_{i}}{e_{i,\mathrm{ref}}}-1$, between the eigenvalue spectrum
$\\{e_{i}\\}_{i=1}^{d}$ of each type of precision matrix estimate and that of
the reference precision matrix, $\\{e_{i,\mathrm{ref}}\\}_{i=1}^{d}$, all
computed using all $n=$2048$$ mock measurements. The dashed grey line
indicates the zero reference. Note that the inverted covariance shrinkage
estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with target
$\mathsfbf{T}^{(2)}$ (solid orange line) is at a constant offset around
$1\text{\,}\mathrm{\char 37\relax}$, as it is very close to the sample
estimate but without the Hartlap factor, which is about
$1\text{\,}\mathrm{\char 37\relax}$ below unity.

Overall, the agreement of the eigenvalue spectra of the estimates with that of
the reference precision matrix $\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{ref}}$
improves as the sample size $n$ increases, with almost perfect agreement when
$n=$2048$$ as all estimates become close to the reference precision matrix
itself. There are several other observations of interest:

  * •

The sample estimates $\mathsfbf{\Pi}$ tend to lower the smallest eigenvalues,
potentially leading to overestimated uncertainties;

  * •

The covariance shrinkage estimates, $\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$,
with the target choice $\mathsfbf{T}^{(2)}$ tend to enlarge the largest
eigenvalues, potentially leading to underestimated uncertainties. However,
this does not occur for the target choice $\mathsfbf{T}^{(1)}$;

  * •

The NERCOME estimates, $\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$, tend to
decrease the largest eigenvalues, potentially leading to overestimated
uncertainties;

  * •

The precision shrinkage estimates $\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}$
tend to slightly decrease the eigenvalues over the entire spectrum with the
exception of the largest eigenvalues, which are more scattered. This could
potentially lead to overestimated uncertainties.

It is also observed that when we use $\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}$ as
target, one out of the $85$ precision shrinkage estimates with $n=$24$$ has
negative eigenvalues, indicated by the grey line dropping towards zero in the
bottom right plot in Fig. 4. We have noted that this precision shrinkage
estimate with negative eigenvalues also has a negative $\hat{\alpha}^{*}$
value. We also note that for $n=$2048$$, the eigenvalues of the covariance
shrinkage estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with target
$\mathsfbf{T}^{(2)}$ are at a constant offset around $1\text{\,}\mathrm{\char
37\relax}$ from the reference eigenvalues, as indicated by the solid orange
line in Fig. 6; this is because for this particular estimate, the shrinkage
intensity parameter $\hat{\lambda}^{*}=$0.002$$ is close to zero, so that it
is effectively the sample estimate itself albeit without the Hartlap factor
that is about $1\text{\,}\mathrm{\char 37\relax}$ below unity.

To make a more direct quantitative comparison, we consider the following loss
function for a generic precision matrix estimate $\hat{\mathsfbf{\Psi}}$:

$\mathop{}\\!L{\bigl{(}\hat{\mathsfbf{\Psi}};\mathsfbf{\Pi}_{\mkern
2.0mu\mathrm{ref}}\bigr{)}}=\mathopen{}\mathclose{{}\left\|\hat{\mathsfbf{\Psi}}^{1/2}\mathsfbf{\Pi}_{\mkern
2.0mu\mathrm{ref}}^{-1}\hat{\mathsfbf{\Psi}}^{1/2}-\mathsfbf{I}_{d}}\right\|_{\mathrm{F}}\,,$
(17)

where $\mathsfbf{I}_{d}$ is the $d$-dimensional identity matrix. Over all
subsets of $n$ mock measurements where $n=$24$$ or $30$, we compute the median
loss function and its uncertainty bounds given by the
$16\text{\,}\mathrm{\char 37\relax}$ and $84\text{\,}\mathrm{\char 37\relax}$
quantiles in Table 3, with additionally a single loss function value for
$n=$2048$$ as there is only one precision matrix estimate of each type using
all available mock measurements.

Table 3:  Median of the loss function $\mathop{}\\!L{\bigl{(}\hat{\mathsfbf{\Psi}}\bigr{)}}$ (equation 17) for each type of precision matrix estimate $\hat{\mathsfbf{\Psi}}$ with different sample sizes $n$, with uncertainty bounds given by the $16\text{\,}\mathrm{\char 37\relax}$ and $84\text{\,}\mathrm{\char 37\relax}$ quantiles. Note that there is only a single loss function value when $n=$2048$$, as there is only one estimate of each type computed from all available mock measurements.  Precision estimation type | $n=$24$$ | $n=$30$$ |  $n=$2048$$
---|---|---|---
Sample | ${6.8}\begin{subarray}{c}+{7.4}\\\ -{2.2}\end{subarray}$ | ${5.4}\begin{subarray}{c}+{1.8}\\\ -{1.2}\end{subarray}$ | –
Covariance shrinkage, target $\mathsfbf{T}^{(1)}$ (empirical) | ${2.4}\begin{subarray}{c}+{1.0}\\\ -{0.3}\end{subarray}$ | ${2.3}\begin{subarray}{c}+{0.6}\\\ -{0.3}\end{subarray}$ |  $0.051$
Covariance shrinkage, target $\mathsfbf{T}^{(2)}$ (analytical) | ${38.3}\begin{subarray}{c}+{25.2}\\\ -{10.5}\end{subarray}$ | ${18.7}\begin{subarray}{c}+{5.3}\\\ -{4.3}\end{subarray}$ |  $0.048$
NERCOME | ${3.6}\begin{subarray}{c}+{1.2}\\\ -{0.5}\end{subarray}$ | ${3.3}\begin{subarray}{c}+{0.5}\\\ -{0.4}\end{subarray}$ |  $0.049$
Precision shrinkage, target $\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(1)}$ (empirical) | ${7.2}\begin{subarray}{c}+{7.8}\\\ -{2.3}\end{subarray}$ | ${5.8}\begin{subarray}{c}+{2.2}\\\ -{1.6}\end{subarray}$ |  $0.064$
Precision shrinkage, target $\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}$ (analytical) | ${6.6}\begin{subarray}{c}+{8.5}\\\ -{3.0}\end{subarray}$ | ${4.7}\begin{subarray}{c}+{1.7}\\\ -{1.1}\end{subarray}$ |  $0.013$

When the sample size is very close to the data vector dimension $d=$18$$, e.g.
$n=$24$$, the covariance shrinkage estimate
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with the empirical target
$\mathsfbf{T}^{(1)}$ gives the smallest loss closely followed by the NERCOME
estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$. The precision shrinkage
estimate $\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}$ with the analytical target
$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}$ has a loss function close to that of
the sample estimate $\mathsfbf{\Pi}$ while the other estimates result in
larger losses. When $n=$30$$ mock measurements are used, most shrinkage
estimates have a lower loss function than the sample estimate, except the
covariance shrinkage estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with the
analytical target $\mathsfbf{T}^{(2)}$ and the precision shrinkage estimate
$\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}$ with the empirical target
$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(1)}$. When $n=$2048$$, all loss function
values are small. Overall, both the covariance shrinkage estimate
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with the empirical target
$\mathsfbf{T}^{(1)}$ and the NERCOME estimate
$\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$ give consistently small loss values
irrespective of the sample size, whereas the shrinkage estimates with an
analytical target have larger loss function values but improve rapidly with
the sample size $n$.

###  4.2 Parameter inference

Since a major objective of precision matrix estimation is to obtain accurate
cosmological parameter constraints in a likelihood analysis, in this section
we consider the log-likelihood function (up to an additive constant) for each
type of precision matrix estimate $\hat{\mathsfbf{\Psi}}$:

$\ln\mathop{}\\!\mathcal{L}{\bigl{(}b,f;\hat{\mathsfbf{\Psi}}\bigr{)}}=\\\
-\frac{1}{2}{\mathopen{}\mathclose{{}\left[\textbf{{P}}_{\mathrm{data}}-\textbf{{P}}_{\mathrm{model}}(b,f)}\right]}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\hat{\mathsfbf{\Psi}}\mathopen{}\mathclose{{}\left[\textbf{{P}}_{\mathrm{data}}-\textbf{{P}}_{\mathrm{model}}(b,f)}\right]\,,$
(18)

where $\textbf{{P}}_{\mathrm{data}}$ is the data vector of power spectrum
monopole and quadrupole, and $\textbf{{P}}_{\mathrm{model}}$ is the model
vector of power spectrum monopole, quadrupole and hexadecapole transformed by
the wide-angle and survey window matrices (see section 3 for details). Note
that for simplicity, we only consider varying the linear bias and growth-rate
parameters $b$ and $f$ with all other parameters fixed to the Planck 2015
cosmology. To perform Bayesian parameter inference, we assume uniform priors
$b\sim\mathcal{U}($0.5$,$3.5$)$ and $f\sim\mathcal{U}($0$,$2$)$, and sample
the posterior distribution using preconditioned Monte Carlo (PMC) implemented
by the pocoMC sampler (Karamanis et al., 2022b, a).

For each type of precision matrix estimate computed from $n=$24$$ mock
measurements, we have obtained $85$ chains of the same number of posterior
samples, which are then combined into a single chain of $170\,000$ parameter
points for each estimate type; for precision matrix estimates from $n=$30$$
mock measurements, $68$ posterior sample chains of the same length have been
generated and combined into a single chain of $136\,000$ parameter points for
each estimate type; and for $n=$2048$$ where all mock measurements are used
for precision matrix estimation, a single chain of at least $2000$ posterior
samples has been generated for each estimate type.

In Figs. 7, 8 and 9, we show the posterior contours for the
$68\text{\,}\mathrm{\char 37\relax}$ and $95\text{\,}\mathrm{\char 37\relax}$
credible regions from each type of precision matrix estimate for sample size
$n=$24$$, $30$ and $2048$ respectively, compared to the reference results
obtained using the reference precision matrix $\mathsfbf{\Pi}_{\mkern
2.0mu\mathrm{ref}}$.

Figure 7:  Posterior constraints in $68\text{\,}\mathrm{\char 37\relax}$ and
$95\text{\,}\mathrm{\char 37\relax}$ credible regions on the linear bias $b$
and growth rate $f$ obtained using different precision matrix estimates
computed from subsets of $n=$24$$ mock measurements (red contours), compared
to the results obtained from the reference precision matrix
$\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{ref}}$ (blue contours). The non-reference
posterior mean parameter estimates are annotated above each subplot with
uncertainties given by the $68\text{\,}\mathrm{\char 37\relax}$ credible
interval.

Figure 8:  The same as Fig. 7 but with precision matrix estimates computed
from subsets of $n=$30$$ mock measurements using the different estimators
discussed in the text.

Figure 9:  The same as Figs. 7 and 8 but with precision matrix estimates
computed from all $n=$2048$$ mock measurements using the different shrinkage
estimators discussed in the text.

In each of these figures, we have also provided the posterior mean estimates
for parameters $b$ and $f$ with the uncertainty bounds given by the
$68\text{\,}\mathrm{\char 37\relax}$ credible interval, compared with the
reference constraints $b=$1.945(18)$$ and $f=$0.617(43)$$. In principle, the
uncertainty in the precision matrix estimate will increase the uncertainties
of the inferred parameters; for the sample precision matrix estimate, it is
possible to derive a multiplicative factor that accounts for this in the
uncertainties of the maximum likelihood estimator (Dodelson & Schneider,
2013), which remains a good approximation for the posterior mean estimator
when the posterior distribution is close to being Gaussian (Percival et al.,
2014). However, it is difficult to derive a similar factor for a generic
shrinkage estimator as its probability distribution may not have an analytical
form. Therefore, instead of applying this correction factor to the inferred
parameter uncertainties, we rely on the fact we have combined posterior sample
chains with multiple precision matrix estimates generated from equivalent
subsets of $n$ mock measurements when $n=$24$$ or $30$, which should reflect
the scatter in the precision matrix estimation.

When the sample size is $n=$24$$, which is small and close to the data vector
dimension $d=$18$$, we see in Fig. 7 that the noisy sample estimate
$\mathsfbf{\Pi}$ significantly worsens the parameter constraints. The
covariance shrinkage estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with the
empirical target $\mathsfbf{T}^{(1)}$ provides constraints much closer to the
reference results, albeit very slightly over-tightened, whereas the analytical
target $\mathsfbf{T}^{(2)}$ leads to overestimated parameter uncertainties and
a few outliers out of the $85$ chains sampled. This explains the protuberance
in the contour boundaries for the covariance shrinkage estimate
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with target choice $\mathsfbf{T}^{(2)}$
in Fig. 7. The NERCOME estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$
similarly gives looser parameter constraints, and in particular there is a
shift in the mean of the posterior distribution of the linear bias $b$,
quantified at the $0.7\text{\,}\mathrm{\text{$\sigma$}}$ level with respect to
the reference results. The precision shrinkage estimate
$\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}$ leads to broader contours that have
a different orientation to the rest, and this gives overestimated
uncertainties on $b$, as well as slightly underestimated uncertainties on $f$
for the analytical target $\mathsfbf{T}^{(2)}$.

With the sample size down to $n=$30$$, improvements are seen across all types
of precision matrix estimates. Again the covariance shrinkage estimate
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with the empirical target
$\mathsfbf{T}^{(1)}$ provides posteriors closest to the reference results,
albeit with slightly over-tightened parameter constraints. The sample estimate
$\mathsfbf{\Pi}$ has the widest parameter bounds, followed by the NERCOME
estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$ and then the other shrinkage
estimates. The shift in the posterior distribution of the linear bias $b$
produced by the NERCOME estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$
remains, although the posterior mean is biased with respect to the reference
results at a reduced level of $0.5\text{\,}\mathrm{\text{$\sigma$}}$. Compared
with sample size $n=$24$$, the difference in the orientation of the posterior
contours between the precision shrinkage estimates $\mathsfbf{\Pi}_{\mkern
2.0mu\mathrm{LS}}$ and the rest is lessened.

When computed from all $2048$ mock measurements, all types of precision matrix
estimates yield posterior constraints in close approximation to those obtained
from using the reference precision matrix $\mathsfbf{\Pi}_{\mkern
2.0mu\mathrm{ref}}$, as seen in Fig. 9. This indicates that given a large
sample size, all estimation methods behave consistently as expected.

To quantify the differences in the posterior constraints with a varying sample
size $n$ (as shown in Figs. 7, 8 and 9), we compute the mean $\hat{b}$ or
$\hat{f}$ and the corresponding standard deviation $\hat{\sigma}$ of the
marginal posterior distributions from the combined chain for each type of
precision matrix estimate; these are then compared against the reference
posterior mean $b_{\mathrm{ref}}$ or $f_{\mathrm{ref}}$ and standard deviation
$\sigma_{\mathrm{ref}}$ extracted from the chain obtained with the reference
precision matrix, as shown in Fig. 10.

Figure 10:  Relative deviations in the posterior mean and standard deviation
of the linear bias $b$ (left panel) and growth rate $f$ (right panel) obtained
from the various precision matrix estimates with a varying sample size $n$,
compared to the reference results $b_{\mathrm{ref}}$ or $f_{\mathrm{ref}}$ and
$\sigma_{\mathrm{ref}}$ obtained using the reference precision matrix. The
dashed grey lines indicate the zero reference.

It is visually clear from Fig. 10 that when $n$ is small, the NERCOME estimate
leads to a noticeable shift in the posterior mean for the linear bias $b$ and
the sample estimate leads to the most overestimated parameter uncertainties,
while both the covariance shrinkage and precision shrinkage estimates produce
results that are closer to the reference.

##  5 Conclusion

The determination of the covariance and precision matrices is vital to
parameter inference in the era of precision cosmology. In this work, we have
focused on a class of estimation methods known as shrinkage, which are used to
obtain the following precision matrix estimates:

  1. 1.

the covariance shrinkage estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$
obtained after inversion with two target choices $\mathsfbf{T}^{(1)}$ and
$\mathsfbf{T}^{(2)}$, where the former is empirically based on the diagonal of
the sample covariance matrix estimate and the latter is analytical (see
section 3.2);

  2. 2.

the NERCOME (nonlinear shrinkage) estimate
$\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$ obtained after inversion (see section
3.3);

  3. 3.

the precision shrinkage estimate $\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}$
with two target choices $\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(1)}$ and
$\mathsfbf{\Pi}_{\mkern 2.0mu0}^{(2)}$, where the former is empirically based
on the eigenvalues of the inverted sample covariance matrix estimate and the
latter is analytical (see section 3.4).

All of these shrinkage estimates, as well as the sample estimate
$\mathsfbf{\Pi}$, are computed for a varying sample size $n$ and compared to
the reference precision matrix $\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{ref}}$
obtained for a large value of $n$.

We have applied these estimation methods to the power spectrum analysis of the
BOSS DR12 data set with data vector dimension $d=$18$$, in conjunction with
$2048$ Patchy mock catalogues for precision matrix estimation. The performance
of the different precision matrix estimates for sample size $n=$24$$, $30$ or
$2048$ is compared against the reference precision matrix constructed using
all mock measurements. We have considered diagnostic metrics such as the
matrix loss function and the eigenvalue spectrum; in addition, we have also
performed Bayesian parameter inference of the linear bias and growth-rate
parameters $b$ and $f$. It has been found that:

  * •

When the sample size is small, the covariance shrinkage estimate
$\hat{\mathsfbf{\Psi}}_{\mathrm{LS}}$ with the empirical target
$\mathsfbf{T}^{(1)}$ provides the lowest loss function values, the closest
eigenvalue spectrum and the best agreement in the posterior parameter
constraints compared to the reference results. With the analytical target
$\mathsfbf{T}^{(2)}$, however, it has a higher loss function value and weaker
parameter constraints.

  * •

The NERCOME estimate $\hat{\mathsfbf{\Psi}}_{\mathrm{NLS}}$ tends to
underestimate the largest precision matrix eigenvalues leading to weaker
parameter constraints, which has also been demonstrated by Gouyou Beauchamps
et al. (2023). It also produces a shift in the posterior mean of the bias
parameter by as much as $0.7\text{\,}\mathrm{\text{$\sigma$}}$ with respect to
the reference results when $n=$24$$.

  * •

The precision shrinkage estimate $\mathsfbf{\Pi}_{\mkern 2.0mu\mathrm{LS}}$
produces inflated parameter constraints with a slightly rotated parameter
degeneracy direction for small sample size $n=$24$$ or $30$, but still
outperforms the sample estimate $\mathsfbf{\Pi}$.

  * •

As expected, when the sample size $n$ is very large, all precision estimates
behave similarly and give consistent posterior constraints.

Based on these configurations and results, we suggest the use of covariance
shrinkage estimation of the precision matrix obtained after inversion. This
serves as a useful guide for the cosmological analyses of forthcoming data
releases of Stage IV galaxy surveys such as DESI and Euclid, where large data
sets are best exploited with more performant precision matrix estimates.

Finally, we note there are many other possible choices of the shrinkage target
not exhausted in this work, and their performance will depend on the precise
configuration of individual cosmological analyses. In particular, although the
precision shrinkage estimate has been found to give weaker parameter
constraints with a rotated parameter degeneracy direction, it needs not be the
case with other unexplored target choices or in a different analysis set-up.
It is expected that with a target matrix that is more structurally similar to
the true precision matrix $\mathsfbf{\Psi}$ (instead of the diagonal targets
considered in section 3.4), e.g. an analytical target with non-linear and
window function effects, these precision shrinkage estimates can be further
improved. Moreover, it would offer new insight to derive correction factors
similar to those in Hartlap et al. (2006) and Percival et al. (2014) to
account for the stochasticity in shrinkage estimates of the precision matrix,
even when they are already much less noisy than the standard sample estimates.
However, this is generally difficult since shrinkage estimates are composite
quantities for which the probability distribution may be analytically
intractable. We leave these considerations to future work.

## Acknowledgements

The authors would like to thank Minas Karamanis for inspiration, and Uroš
Seljak, Richard Grummit, Benjamin Joachimi, Julian Bautista and Taras Bodnar
for helpful discussions.

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement 853291). MJL acknowledges support from the School of Physics
& Astronomy at the University of Edinburgh. FB acknowledges the support of the
Royal Society through a University Research Fellowship.

This work has made use of publicly available Python packages NumPy (Harris et
al., 2020), SciPy (Virtanen et al., 2020), Matplotlib (Hunter, 2007), GetDist
(Lewis, 2019), nbodykit (Hand et al., 2018) and pocoMC (Karamanis et al.,
2022a).

## Data Availability

The data underlying this article will be shared on request to the
corresponding author. The source data are available at
fbeutler.github.io/hub/deconv_paper, and any generated data, including the
code used to derive further data, are available at
github.com/marnixlooijmans/shrinkage-estimation-paper.

## REFERENCES

  * Abadir et al. (2014) Abadir K. M., Distaso W., Žikeš F., 2014, J. Econom., 181, 165
  * d’Amico et al. (2020) d’Amico G., Gleyzes J., Kokron N., Markovic K., Senatore L., Zhang P., Beutler F., Gil-Marín H., 2020, J. Cosmology Astropart. Phys., 2020, 005
  * BOSS Collaboration et al. (2015) BOSS Collaboration et al., 2015, ApJS, 219, 12
  * Beutler & McDonald (2021) Beutler F., McDonald P., 2021, J. Cosmology Astropart. Phys., 2021, 031
  * Bodnar et al. (2016) Bodnar T., Gupta A. K., Parolya N., 2016, J. Multivar. Anal., 146, 223
  * DESI Collaboration (2023) DESI Collaboration 2023, preprint (arXiv2306.06308)
  * DESI Collaboration (2024) DESI Collaboration 2024, AJ, 167, 62
  * Dodelson & Schneider (2013) Dodelson S., Schneider M. D., 2013, Phys. Rev. D, 88, 063537
  * Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, ApJ, 496, 605
  * Escoffier et al. (2016) Escoffier S. et al., 2016, preprint (arXiv1606.00233)
  * Friedrich et al. (2015) Friedrich O., Seitz S., Eifler T. F., Gruen D., 2015, MNRAS, 456, 2662
  * Gouyou Beauchamps et al. (2023) Gouyou Beauchamps S., Baratta P., Escoffier S., Gillard W., Bel J., Bautista J., Carbone C., 2023, preprint (arXiv1110.3193)
  * Grieb et al. (2016) Grieb J. N., Sánchez A. G., Salazar-Albornoz S., Vecchia C. D., 2016, MNRAS, 457, 1577
  * Hamilton et al. (2006) Hamilton A. J., Rimes C. D., Scoccimarro R., 2006, MNRAS, 371, 1188
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, AJ, 156, 160
  * Harris et al. (2020) Harris C. R. et al., 2020, Nature, 585, 357
  * Hartlap et al. (2006) Hartlap J., Simon P., Schneider P., 2006, A&A, 464, 399
  * Hunter (2007) Hunter J. D., 2007, Comput. Sci. Eng., 9, 90
  * Joachimi (2017) Joachimi B., 2017, MNRAS, 466, L83
  * Kaiser (1987) Kaiser N., 1987, MNRAS, 227, 1
  * Karamanis et al. (2022a) Karamanis M., Nabergoj D., Beutler F., Peacock J. A., Seljak U., 2022a, J. Open Source Softw., 7, 4634
  * Karamanis et al. (2022b) Karamanis M., Beutler F., Peacock J. A., Nabergoj D., Seljak U., 2022b, MNRAS, 516, 1644
  * Kitaura et al. (2016) Kitaura F. S. et al., 2016, MNRAS, 456, 4156
  * Lam (2016) Lam C., 2016, Ann. Stat., 44, 928
  * Laureijs et al. (2011) Laureijs R. et al., 2011, preprint (arXiv1110.3193)
  * Ledoit & Wolf (2003) Ledoit O., Wolf M., 2003, J. Empir. Financ., 10, 603
  * Ledoit & Wolf (2004) Ledoit O., Wolf M., 2004, Journal of Multivariate Analysis, 88, 365
  * Ledoit & Wolf (2012) Ledoit O., Wolf M., 2012, Ann. Stat., 40, 1024
  * Lewis (2019) Lewis A., 2019, preprint (arXiv1910.13970)
  * Li et al. (2019) Li Y., Singh S., Yu B., Feng Y., Seljak U., 2019, J. Cosmology Astropart. Phys., 2019, 016–016
  * Mohammad & Percival (2022) Mohammad F. G., Percival W. J., 2022, MNRAS, 514, 1289
  * O’Connell & Eisenstein (2019) O’Connell R., Eisenstein D. J., 2019, MNRAS, 487, 2701–2717
  * Percival et al. (2014) Percival W. J. et al., 2014, MNRAS, 439, 2531
  * Perko et al. (2016) Perko A., Senatore L., Jennings E., Wechsler R. H., 2016, preprint (arXiv1610.09321)
  * Philcox et al. (2019) Philcox O. H. E., Eisenstein D. J., O’Connell R., Wiegand A., 2019, MNRAS, 491, 3290–3317
  * Planck Collaboration (2016) Planck Collaboration 2016, A&A, 594, A13
  * Pope & Szapudi (2008) Pope A. C., Szapudi I., 2008, MNRAS, 389, 766
  * Reid et al. (2016) Reid B. et al., 2016, MNRAS, 455, 1553
  * Schäfer & Strimmer (2005) Schäfer J., Strimmer K., 2005, Stat. Appl. Genet. Mol. Biol., 4
  * Sellentin & Heavens (2016) Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132
  * Semenaite et al. (2022) Semenaite A. et al., 2022, MNRAS, 512, 5657
  * Sugiyama et al. (2020) Sugiyama N. S., Saito S., Beutler F., Seo H. J., 2020, MNRAS, 497, 1684
  * Taylor & Joachimi (2014) Taylor A., Joachimi B., 2014, MNRAS, 442, 2728
  * Taylor et al. (2013) Taylor A., Joachimi B., Kitching T., 2013, MNRAS, 432, 1928
  * de la Torre et al. (2013) de la Torre S. et al., 2013, A&A, 557, A54
  * Virtanen et al. (2020) Virtanen P. et al., 2020, Nature Methods, 17, 261
  * Wadekar & Scoccimarro (2020) Wadekar D., Scoccimarro R., 2020, Phys. Rev. D, 102, 123517
  * Wadekar et al. (2020) Wadekar D., Ivanov M. M., Scoccimarro R., 2020, Phys. Rev. D, 102, 123521