# IOP4, the Interactive Optical Photo-Polarimetric Python Pipeline

Juan Escudero Pedrosa Instituto de Astrofísica de Andalucía, Glorieta de la
Astronomía, s/n. Granada, 18008, Spain  Iván Agudo Instituto de Astrofísica de
Andalucía, Glorieta de la Astronomía, s/n. Granada, 18008, Spain  Daniel
Morcuende Instituto de Astrofísica de Andalucía, Glorieta de la Astronomía,
s/n. Granada, 18008, Spain  Jorge Otero-Santos Instituto de Astrofísica de
Andalucía, Glorieta de la Astronomía, s/n. Granada, 18008, Spain  Giacomo
Bonnoli INAF Osservatorio Astronomico di Brera, Via E. Bianchi 46. Merate
(LC), 23807, Italy  Vilppu Piirola  Department of Physics and Astronomy,
University of Turku, Vesilinnantie 5. Turku, FI-20014, Finland  César Husillos
Geological and Mining Institute of Spain (IGME-CSIC), Calle Ríos Rosas 23,
E-28003, Madrid, Spain  Instituto de Astrofísica de Andalucía, Glorieta de la
Astronomía, s/n. Granada, 18008, Spain  Mabel Bernardos  Instituto de
Astrofísica de Andalucía, Glorieta de la Astronomía, s/n. Granada, 18008,
Spain  Rubén López-Coto Instituto de Astrofísica de Andalucía, Glorieta de la
Astronomía, s/n. Granada, 18008, Spain  Alfredo Sota  Instituto de Astrofísica
de Andalucía, Glorieta de la Astronomía, s/n. Granada, 18008, Spain  Victor
Casanova  Instituto de Astrofísica de Andalucía, Glorieta de la Astronomía,
s/n. Granada, 18008, Spain  Fran J. Aceituno  Instituto de Astrofísica de
Andalucía, Glorieta de la Astronomía, s/n. Granada, 18008, Spain  Pablo
Santos-Sanz  Instituto de Astrofísica de Andalucía, Glorieta de la Astronomía,
s/n. Granada, 18008, Spain

###### Abstract

IOP4 is a pipeline to perform photometry and polarimetry analysis of optical
data from Calar Alto (CAHA) and Sierra Nevada (OSN) observatories. IOP4
implements Object Relational Mapping (ORM) to seamlessly integrate all
information about the reduction and results in a database which can be used to
query and plot results, flag data and inspect the reduction process in an
integrated fashion with the whole pipeline. It also ships with an already
built-in web interface which can be used out of the box to browse the database
and supervise all pipeline processes. It is built to ease debugging and
inspection of data. Reduction from five different instruments are already
implemented: RoperT90, AndorT90 and DIPOL (at OSN 0.9m telescope), AndorT150
(OSN 1.5m telescope) and CAFOS (CAHA 2.2m telescope). IOP4’s modular design
allows for easy integration of new observatories and instruments, and its
results have already featured in several high-impact refereed publications. In
this paper we describe the implementation and characteristics of IOP4.

Data analysis — Photometry — Polarimetry – Astronomy databases – Open source
software

††facilities: Instituto de Astrofísica de Andalucía (IAA-CSIC), Observatorio
de Sierra Nevada (OSN), Observatorio de Calar Alto (CAHA)††software:  Django
(Django Software Foundation 2023), SQLite (SQLite Developers 2023), Vue.js
(Vue Developers 2023), numpy (Harris et al. 2020), astropy (Astropy
Collaboration et al. 2013, 2018, 2022), photutils (Bradley et al. 2023),
astrometry.net (Lang et al. 2010)

tablenum SIXtablenum

##  1 Introduction

Optical photo-polarimetric observational programs, especially those dedicated
to monitoring, can regularly produce large quantities of data that can take
considerable time and effort to be managed and reduced. The effort necessary
to produce good quality results extends beyond the use of automatic tools and
can include a human-supervised iterative process of debugging the reduction
and comparing the employed methods and results with those of different
programs. The use of automatic tools is therefore necessary, however, in many
instances they obscure the process of reduction and intermediate results,
making the debugging of any problem in the results a hard task.

IOP4 implements Object Relational Mapping using Django’s ORM system. This
allows to transparently keep the database schema up to date with the models
used by the pipeline without any need to mess with the underlying SQL queries.
The choice of Django’s as ORM backend allows to seamlessly use the rest of
Django Framework to serve the results, including but not limited to its admin
interface to inspect the database, and Django’s debug web server. These tools
are all written in Python and packaged for pip and conda, a programming
language and distributions that many astronomers are already familiarized
with. This also makes IOP4 a multi-platform software, compatible both with
macOS and Linux.

IOP4 builds on top of existing technologies, some of them already cited:
Django (ORM and web application framework), SQLite (default database backend),
astrometry.net (blind astrometric calibration, Lang et al. 2010), Bokeh (high
quality and interactive plots in the web interface), Vue.js (single page web
application), Quasar (user interface components) and JS9 (a ds9 web port for
interactive FITS visualization). Many other open-source packages are used, the
complete list can be found in the setup.toml file. All of them are
automatically installed together with IOP4.

Hardware requirements of IOP4 are low, even compared today to modern day
consumer-end laptops. The blind astrometric calibration using the
astrometry.net solver and its index files can take as much as 35GB with the
default configuration, although this requirement can be lowered. SSD storage
is also recommended, since it needs to read significant amounts of data
especially for the astrometric calibration. Although it is able to run on a
single-core, significant speed improvements can be gained from using up to 20
cores, the recommended maximum with the default configuration and database
(DB) backend. Higher levels of concurrency are possible by tuning the DB
configuration. A good internet connection can speed up the first execution of
IOP4, including tests, since it will need to download astrometry index files.

IOP4 facilitates both automatic reduction and manual inspection of the
procedure. The iop4 command provides several options to select different
epochs and files and automatically process them. The iop4lib allows any user
to write its own custom reduction scripts, and provides the tools to invoke
any part of the reduction procedure and inspect and manipulate any object in
the database from an interactive terminal, a Python script or a Jupyter
Notebook. The iop4api implements several API (Application Public Interface)
endpoints and web applications that allow end users to easily query results,
produce plots, flag bad data and inspect any object. The portal can be used
standalone as in the ready for use iop4site project, or integrated into other
sites to be served to the general public.

##  2 General reduction procedure

The automatic reduction procedure generally starts with a simple invocation of
the iop4 command. The main script can be requested to select epochs in the
local archive or to discover and download epochs in the remote telescope
repositories. For each epoch, the script goes through the usual steps of
photo-polarimetric reduction:

  1. 1.

Classification of raw science images: This includes their type (bias, darks,
flats and science images), discovering the instrument and type of observation
(photometry or polarimetry), and translation of standard and non-standard
keywords (exposure time, band, datetime, rotator angles, objects).

  2. 2.

Creation of master calibration frames: Images of each type (bias, darks and
flats) are grouped and merged together to create the master calibration nights
available for each night.

  3. 3.

Reduction and calibration of science images: Raw images are applied the
corresponding master calibration frames. The reduced images are also given a
correct World Coordinate System (WCS) in their header after astrometric
calibration.

  4. 4.

Computation of photo-polarimetric results: At the moment, the procedure
implements relative photometry using known calibrators in the field, and
polarimetry both for half-wave ($\lambda/2$) retarder based polarizers and
polarized filter wheels.

  5. 5.

Post-processing of results: This includes correction of magnitude and degree
of polarization to account for the host contribution (Nilsson et al. 2007),
and the possibly needed transformation to a standard photometric system, if
the instrument used for the observations does not have one.

##  3 Data organization

IOP4 data directory structure follows the typical hierarchical schema shown in
Fig. 1. In this schema, all raw data is stored and isolated under a single
folder (raw/), which allows to set up a local archive of the original data
without any modifications for long-term conservation, and to set up the
necessary permissions to protect and share it with other users independently
of the rest of IOP4-created files. Under the raw directory, data is organized
first by telescope and then by night of observation. Other files such as built
master calibration frames and reduced images are stored separately. Also
auxiliary images such as automatically built previews, finding charts, summary
plots, etc, which are too heavy to be stored in a database, are stored under
different folders.

Figure 1: Directory structure of the IOP4 local archive.

The database schema is shown in Fig. 2. IOP4 implements object-relational
mapping (ORM). In ORM, objects in a object-oriented programming language
(Python classes) are mapped to tables in a relational database, while
instances of these objects correspond to rows in each of the tables. This
enormously simplifies the interaction with the database, removing the need for
writing SQL queries and manipulating the database to keep its structure and
content updated.

The default database backend is SQLite111https://www.sqlite.org/., a in-
process library that implements a self-contained, serverless, zero-
configuration, and transactional SQL database; and uses the Django
ORM222https://docs.djangoproject.com/. to interact with it. Customizations and
changes to the DB schema can be written in the code (e.g, adding an attribute
to the photo-polarimetric result model) and they will be automatically
propagated to the DB schema by the Django migration system.

Figure 2: Database schema of IOP4. ForeignKey relationships relate one or
several instances of model A to a single instance of model B. In ManyToMany
relationships, multiple instances of model A and B are linked together. The
latter does not show a field in the tables of this diagram, since the
relationship is established through a hidden table omitted here. Some fields
of the AstroSource model have been also omitted to save space (corresponding
to literature magnitudes in other bands, e.g. mag_B, etc). An arrow in any
direction signify multiple instances being linked to the origin (e.g. several
RawFit(s) are linked to one Epoch).

##  4 Astrometric calibration

Astronomical images are usually distributed in the FITS format (Pence et al.
2010). The format allows for a World Coordinate System (WCS, Greisen &
Calabretta 2002) to be incorporated in the metadata or header section of the
FITS file. However, most raw science images from telescopes do not include
this precise information and need to be calibrated. For most images, which
have a wide field of view ($>$7\text{\,}\mathrm{\SIUnitSymbolArcminute}$$)
this is done by a local solver333The local solver is integrated as an external
Python module dependency in the pyproject.toml file and is installed next to
IOP4 automatically. The module is a wrapper around the native C functions of
the astrometry.net library. https://github.com/neuromorphicsystems/astrometry
based on the astrometry.net library
solver444https://github.com/dstndstn/astrometry.net.

For observations with Ordinary (O) and Extraordinary (E) images (e.g. CAFOS
and DIPOL imaging polarimetric observations), astrometric calibration involves
a previous step of separating the detected sources in pairs. The source
pairing is done by finding the most common distance between all pairs of
sources in the image, which results in two distributions like those in Fig. 3.
Then, the separation can proceed by looking at which pairs are at the right
distance (Fig. 4). Without any constraint, this process is not completely
error free, specially for images with few detected sources. However, since the
distance between pairs is usually a stable property of the instrument, the
known distance between pairs (obtained by the unconstrained pairing) for an
instrument can be used as initial input, which improves the success rate up to
$100\text{\,}\mathrm{\char 37\relax}$.

Figure 3: Distribution of the distances between all possible pairs in a
imaging polarimetric frame. The two peaks correspond to the distance between
the ordinary and extraordinary sources in the image. The image and the
corresponding found pairs can be found in Fig. 4. The orange and blue
distributions correspond to the distances in X and Y axis respectively. The
image corresponds to a CAFOS photo-polarimetric observation of the BL Lacertae
field. Figure 4: Paired sources in an example of CAFOS imaging polarimetric
observation of the BL Lacertae field. The distance used for pairing was
obtained from Fig. 3. The resulting calibrated image can be found in Fig. 5.
Ordinary and Extraordinary sources are circled with the same color.

The detected sources in the image, or the ordinary set of pairs for images
with pairs, can then be used as input for the local astrometry solver (Lang et
al. 2010). The solver compares invariant hashes computed from the detected
source positions to pre-computed hashes from astronomical catalogs. Although
it does not require an initial guess of the position nor the pixel size (blind
solving), the speed and accuracy of the solution is greatly improved by
providing the known pixel size for the instrument and a hint position,
obtained from the header of the images. It returns a list of matches and their
corresponding log-odds. Its success is dependent on the source detection step,
and therefore several attempts are made with different detection parameters
(such as signal threshold) until a good match is found. The resulting WCS is
written to the header of the reduced FITS file. The WCS for the second or
extraordinary set of pairs is directly built from the first one by
translation, and written next to it in the same header. An example of the
result can be found in Fig. 5, which corresponds to the pairs in Figs. 3 and
4.

Figure 5: Calibrated CAFOS photo-polarimetric image of the BL Lacertae field.
The positions of the ordinary image of BL Lacertae (labeled as 2200+420) and
its calibrators (B, C, H) are indicated with a circle. The position of the
extraordinary images are represented with an ’x’.

###  4.1 A quad hash for DIPOL polarimetry images

To reduce the disk capacity requirements of DIPOL polarimetry images, only a
subframe of the full field of DIPOL camera is saved (Sect. 6.3). The reduced
field of view does not usually contain enough stars to perform the astrometric
calibration using the default solver. In fact, in many cases, only the O and E
images of the target source are visible in the image. For the case of images
in which only 1-2 pairs of sources appear, the astrometric calibration can be
as simple as using the central pair as a reference for the creation of the WCS
with the known angle of the RA-DEC grid. However, this process is error-prone
when more than two pairs of sources appear in the image, as the likelihood of
choosing the wrong pair increases. Even for the case of high number of pairs,
the size of the subframe
($$2.5\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$2.0\text{\,}\mathrm{\SIUnitSymbolArcminute}$$)
is too small even for the smallest skymarks of the default astrometry.net
index files555http://astrometry.net/doc/readme.html. The problem of choosing
the right pair of sources as the reference can be solved by comparing the
subframe of the polarimetry observation with the central part of a calibrated
photometry (full frame) field. To this end, we have implemented a hash
algorithm loosely based on Lang et al. (2010). The proposed hash is invariant
under rotation and reflection, although not under scaling. An example of this
algorithm at work can be found in Fig. 6. It takes the ten brightest sources
in each image (including E and O images) and compares all the quads between
them. The best matching quads are used for identifying the target star.

Figure 6: DIPOL polarimetry (left) and photometry (right) frames of BL
Lacertae observations. Only the shown subframe
($$2.5\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$2.0\text{\,}\mathrm{\SIUnitSymbolArcminute}$$)
is usually saved to disk during polarimetric observations with DIPOL, which
prevents blind solving of the image. The full photometry field
($$9.2\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$6.3\text{\,}\mathrm{\SIUnitSymbolArcminute}$$)
is saved, however, and the correspondence between both images can be found by
comparing quads of sources through a hashing algorithm. This allows us to
automatically distinguish the O and E images of our target source in the
field.

##  5 High-level web interface

As one of its main goals, IOP4 also provides an API to interact with the data
and a web interface to act as a client. Both are provided by the iop4api
Django application. It defines API end points to query results, produce
interactive web-based plots and flag data. It also integrates the iop4admin
site, a customized Django-admin. The admin site allows us to inspect any
models in the database (Fig. 2).

The IOP4 catalog is part of the database. IOP4 uses a single, unified catalog
that gathers all information about the sources of interest, including the
target sources, calibrators, relationships between them, comments for
observers, etc. Editing of the catalog can be done through IOP4 as with any
other model, or through the web interface, and the changes take immediate
effect on the reduction process of the pipeline.

Figure 7: Different ways of interacting with IOP4. The pipeline (through the
iop4 command, Python script or Jupyter notebooks) can directly interact with
the data aided by iop4lib. The iop4api application exposes a series of API
endpoints that allow the client to authenticate itself, query, plot and flag
data, and explore catalog and logs from a web browser through a Single-Page
Application (SPA). The iop4admin site independently provides a way to interact
and edit all models in the database.

The situation is illustrated in Fig. 7. The exposed API endpoints are used by
the single-page web application (SPA) to authenticate the client, query, plot
and flag data, and explore logs and catalog. The SPA provides a viewer for the
colored logs and allows filtering by logging level (debug, info, warning, and
error) and string searching. The interactive plot is built and serialized in
the server using the Bokeh666https://docs.bokeh.org/en/latest/ (Bokeh
Development Team 2023) Python library, sent to the client and rendered by the
BokehJS library. The interactive plot can be used to directly flag the data.
TabulatorJS777https://tabulator.info/ is used to display the results, and
allows filtering, selecting columns and exporting to several data formats such
as CSV. The SPA itself is built in html, css and javascript using Vue.js
framework as a standalone script to avoid the build step and allow IOP4 to be
installed and used by the astrophysics community in a familiar way (through
pip or conda).

The web application provides an auxiliary tool to facilitate the addition to
the catalog of calibrators for sources with no known previously documented
calibrators. A search for standard stars with constant brightness is performed
within the PanSTARRS888https://catalogs.mast.stsci.edu/panstarrs/ catalog. For
this, we filter stars within the FoV of all instruments that have a relatively
large amount of observations available (typically $N\geq 10$), with a standard
deviation of their aperture SDSS $\mathrm{gri}$ magnitudes $<0.01$. In order
to use as calibrators stars that do not have the risk of saturating the
images, we also restrict the search to targets with magnitudes between $13$
and $18$, typically. With these considerations, we retrieve the $\mathrm{gri}$
aperture magnitudes of non-variable stars that will be used as calibrators.
Due to the different photometric system filter used by the PanSTARRS database
(based on $grizy$ filters) and that from the instruments implemented in IOP4
(generally equipped with standard Johnson-Cousins filters) a conversion
between photometric systems is needed. Three transformations are currently
implemented in IOP4: Jordi et al. (2006), Jester et al. (2005) and Lupton
(2005). As explained by these authors, these transformations are suitable for
stars, with the caveat of Jester et al. (2005) being only suitable for stars
with $R_{C}-I_{C}<1.15$. All three transformations have been found to be
compatible for the calibrators added following this procedure. IOP4 implements
by default the transformations from Lupton (2005).

##  6 Current instruments

Five instruments from three different telescopes are already implemented in
IOP4: RoperT90, AndorT90 and DIPOL (at OSN 0.9m telescope), AndorT150 (OSN
1.5m telescope) and CAFOS (CAHA 2.2m telescope). The modular design of IOP4
allows easily integrating new instruments and observatories. The Telescope
base class provides the skeleton over which to implement new telescopes. It
implements the necessary methods to query, download (e.g. through the
FTPArchiveMixin class) and perform the initial classification of observing
epochs.

Integrating a new instrument is as easy as subclassing the Instrument
subclass. The new subclass must provide the necessary information to identify
the instrument, and possibly implement methods to translate non-standard
keywords (if any), extract position and size hints for astrometry or perform a
custom reduction process.

###  6.1 RoperT90, AndorT90, AndorT150

IOP4 implements data reduction from the current CCD (Andor ikon-L) cameras at
the 0.9m and 1.5m telescopes in Sierra Nevada Observatory (OSN). It also
implements the old Roper (VersArray) cameras (which were installed until
October 2021 and July 2018 respectively). The AndorT90 instrument mounted at
the Nasmyth east focus of the T090 telescope has a
$$13.2\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$13.2\text{\,}\mathrm{\SIUnitSymbolArcminute}$$
field of view and pixel size of
$0.387\text{\,}\mathrm{\SIUnitSymbolArcsecond}{px}^{-1}$. The AndorT150
instrument at the Nasmyth west focus of the T150 has similar characteristics,
with a
$$7.92\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$7.92\text{\,}\mathrm{\SIUnitSymbolArcminute}$$
field of view and pixel size of
$0.232\text{\,}\mathrm{\SIUnitSymbolArcsecond}{px}^{-1}$. They are cooled down
down to $-80\text{\,}\mathrm{\SIUnitSymbolCelsius}$ without need of liquid
nitrogen, and can be further cooled down to
$-100\text{\,}\mathrm{\SIUnitSymbolCelsius}$ using liquid refrigerant. The low
temperatures make the use of dark current calibration frames unnecessary.
Apart from the usual band filter wheels, a polarized filter wheel is available
at the T090 and the T150 that allows polarimetry measurements by taking series
of four images at varying polarized angles in $45\text{\,}\deg$ steps. From
these, the total flux and raw Stokes parameters are computed as

$F=\langle f\rangle=\frac{1}{4}\sum_{i}{f_{i}}$ (1)

and

$\displaystyle q_{\text{raw}}$ $\displaystyle=$
$\displaystyle\frac{f_{0}-f_{90}}{f_{0}+f_{90}}$ (2) $\displaystyle
u_{\text{raw}}$ $\displaystyle=$
$\displaystyle\frac{f_{45}-f_{-45}}{f_{45}+f_{-45}}~{},$ (3)

where $f_{i}$ are the fluxes at each rotator angle. The instrumental
polarization is corrected by applying an offset

$\displaystyle q_{c}$ $\displaystyle=$ $\displaystyle
q_{\text{raw}}-q_{\text{inst}}$ (4) $\displaystyle u_{c}$ $\displaystyle=$
$\displaystyle u_{\text{raw}}-u_{\text{inst}}$ (5)

and a rotation

$\displaystyle q$ $\displaystyle=$ $\displaystyle
q_{c}\cos{\left(2\chi_{\text{inst}}\right)}-u_{c}\sin{\left(2\chi_{\text{inst}}\right)}$
(6) $\displaystyle u$ $\displaystyle=$ $\displaystyle
u_{c}\sin{\left(2\chi_{\text{inst}}\right)}+u_{c}\cos{\left(2\chi_{\text{inst}}\right)}~{},$
(7)

from which the linear polarization degree and polarization angle are computed
as

$\displaystyle p$ $\displaystyle=$ $\displaystyle\sqrt{q^{2}+u^{2}}$ (8)
$\displaystyle\chi$ $\displaystyle=$
$\displaystyle\frac{1}{2}\arctan{\left(u,q\right)}~{}.$ (9)

###  6.2 CAFOS

The CAFOS (Calar Alto Faint Object Spectrograph) instrument, mounted on the
2.2m telescope at Calar Alto Observatory, provides imaging polarimetry
capabilities. It is equipped with a Wollaston prism and a rotable $\lambda/2$
retarder plate that provides two polarized images separated
$18\text{\,}\mathrm{\SIUnitSymbolArcsecond}$. Only a $$800$\times$800$$
subframe of the full $$2048$\times$2048$$ CCD chip is typically used, with a
field of view of
$$34\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$34\text{\,}\mathrm{\SIUnitSymbolArcminute}$$
and a pixel size of $0.530\text{\,}\mathrm{\SIUnitSymbolArcsecond}{px}^{-1}$.
The CCD chip is cooled to temperatures lower than
$-100\text{\,}\mathrm{\SIUnitSymbolCelsius}$, making the dark current
negligible. Several filters are available, while polarimetry observations are
usually taken in Johnson R. The ordinary (O) and extraordinary (E) images at
each angle are used to compute the total flux and the reduced Stokes
parameters as (Zapatero Osorio et al. 2005)

$\displaystyle F=\frac{1}{N}\sum_{i}\frac{f_{E,i}+f_{O,i}}{2}$ (10)

and

$\displaystyle R_{Q}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{f_{O,0}/f_{E,0}}{f_{O,45}/f_{E,45}}}$ (11)
$\displaystyle R_{U}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{f_{O,22}/f_{E,22}}{f_{O,67}/f_{E,67}}}~{},$ (12)

where $f_{O,i}$ and $f_{E,i}$ are the fluxes of the ordinary and extraordinary
images of the source at each angle $i$. From these, the Stokes parameters are
reconstructed as

$\displaystyle Q_{I}$ $\displaystyle=$ $\displaystyle\frac{R_{Q}-1}{R_{Q}+1}$
(13) $\displaystyle U_{I}$ $\displaystyle=$
$\displaystyle\frac{R_{U}-1}{R_{U}+1}~{},$ (14)

and the polarization degree and polarization angles from Eqs. 8 and 9. While
instrumental polarization is negligible for CAFOS, the polarization angle
still needs to be corrected according to

$\chi=\chi-CPA~{},$ (15)

where CPA is the Zero Polarization Angle.

###  6.3 DIPOL

The DIPOL-1 polarimeter is thoroughly described in Piirola, V. et al. (2020)
and Otero-Santos et al. (2024). It is based on a $\lambda/2$ retarder plate
attached to a rotator and a high readout speed CMOS camera. Installed at the
OSN-T090 the instrument has a field of view of
$$9.2\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$6.3\text{\,}\mathrm{\SIUnitSymbolArcminute}$$
and a pixel scale of $0.134\text{\,}\mathrm{\SIUnitSymbolArcsecond}{px}^{-1}$.
Cycles of 16 images are typically taken with varying rotator angles at
$22.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ steps. For polarimetry
observations, usually only a subframe of size
$$2.5\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$2.0\text{\,}\mathrm{\SIUnitSymbolArcminute}$$
is saved, greatly reducing used disk space. Dark current calibration frames
are necessary. High throughput sharp cutoff R, G and B filters are available.
The computation of the polarimetric results is done following Patat &
Romaniello (2006). The Stokes parameters are computed as

$\displaystyle Q_{\text{raw}}$ $\displaystyle=$
$\displaystyle\frac{2}{N}\sum_{i}F_{i}\cos{\left(\frac{\pi}{2}i\right)}$ (16)
$\displaystyle U_{\text{raw}}$ $\displaystyle=$
$\displaystyle\frac{2}{N}\sum_{i}F_{i}\sin{\left(\frac{\pi}{2}i\right)}~{},$
(17)

where the coefficients $F_{i}$ are computed as

$F_{i}=(F_{O,i}-F_{E,i})/(F_{O,i}+F_{E,i})~{},$ (18)

$F_{O,i}$ and $F_{E,i}$ being the fluxes of the ordinary and extraordinary
images of the source for each rotator angle $i$. The instrumental polarization
is corrected with an offset

$\displaystyle Q$ $\displaystyle=$ $\displaystyle
Q_{\text{raw}}-Q_{\text{inst}}$ (19) $\displaystyle U$ $\displaystyle=$
$\displaystyle U_{\text{raw}}-U_{\text{inst}}~{}.$ (20)

Then, $p$ and $\chi$ are obtained per Eqs. 8, 9 and 15. Photometric (full
frame) observations are also performed. The conversion of the magnitudes in
the Baader $\mathrm{RGB}$ filters to standard Johnson-Cousins
$\mathrm{UBVR_{c}I_{c}}$ will be treated in a separate paper (in preparation).

##  7 Development and CI

IOP4 is open-source, its code hosted at https://github.com/juanep97/iop4. The
repository contains all the necessary code to run IOP4 both as a program (iop4
script) and as a library (iop4lib), plus the iop4api Django application and
the customized iop4admin site. Both can be readily used through the Django
debug server and the provided iop4site, or they can be integrated into another
site and deployed for public use (see Serving IOP4 in production in IOP4
documentation). To ease the development process and the maintainability of the
code, we implemented Continuous Integration and Continuous Deployment (CI/CD)
workflows to perform automatic testing, build the documentation and facilitate
the delivery of IOP4. The source code includes a test suit, the test dataset
is freely available at https://vhega.iaa.es/iop4/. The size of the
astrometry.net solver index files makes inviable using GitHub-hosted ordinary
runners for CI. Instead, a self-hosted runner provided by the
VHEGA999https://vhega.iaa.es/ group with self-provisioning capabilities using
garm101010https://github.com/cloudbase/garm automatically runs tests on pull
requests (PR) and merge commits to the main branch on isolated containers that
already provide the test dataset and astrometry index files. Anyone can run
the test locally in their computers using
pytest111111https://docs.pytest.org/. The documentation can be built using
Sphinx121212https://www.sphinx-doc.org/, and contains several notebook
examples. MyST-NB131313https://myst-nb.readthedocs.io/ and
Jupytext141414https://jupytext.readthedocs.io/ allow using the percent format
for the notebooks, an human-readable format easily integrable with version
control software. The example notebooks are automatically run when building
the documentation and use the test dataset. The documentation is also
automatically built and deployed to GitHub Pages as part of the CI. New
releases of IOP4 are automatically deployed to the public PyPi software
repository151515https://pypi.org/project/iop4/.

##  8 Conclusions

IOP4 is an open source, interactive photo-polarimetric pipeline written in
Python. Its results have already featured in a number of refereed and high
impact publications (e.g. Middei et al. (2023a); Peirson et al. (2023); Middei
et al. (2023b); Ehlert et al. (2023); Di Gesu et al. (2023); Marshall et al.
(2023); Kim et al. (2024); Otero-Santos et al. (2024); among others),
following the footsteps of its predecessor IOP3161616Husillos et al. (2021)
(e.g. Liodakis et al. (2022), etc), and has provided data to many on-going
studies. Most of this data comes from the MAPCAT (Agudo et al. 2012) and TOP-
MAPCAT programs, having reduced more than 500GB of data from the these
programs as of February 2024. The former was running at Calar Alto from 2007
to 2018, the latter has been running from 2018 to the present day both at OSN
and CAHA. Both programs are focused on the photo-polarimetric monitoring of
blazars, combining regular observations with targets of opportunity (ToO).
Aided by its parallel processing capabilities, it routinely downloads, reduces
and serves results from a full night of observation of these programs in less
than half an hour. Its speed makes it suitable to be used as a real-time
analysis tool. Moreover, IOP4 has contributed to the publication of several
Astronomer Telegrams (such as Otero-Santos et al. 2023a, b) thanks to the
promptness of its results and the ease of use.

Development of IOP4 is on-going. Future releases of IOP4 might include several
other instruments and reduction methods. As of version v1.1.2, IOP4 provides a
prefabricated night summary script that sends the results of observations to
subscribed users. Future versions might include a more advanced alert and
trigger system. In any case, the iop4lib already allows the users to create
their own scripts, for alerts or any other purpose.

The project welcomes any interested user to participate in IOP4 development
and request or contribute to the implementation of new instruments and
features.

The IAA-CSIC team acknowledges financial support from the Spanish ”Ministerio
de Ciencia e Innovación” (MCIN/AEI/ 10.13039/501100011033) through the Center
of Excellence Severo Ochoa award for the Instituto de Astrofísica de
Andalucía-CSIC (CEX2021-001131-S), and through grants PID2019-107847RB-C44 and
PID2022-139117NB-C44. P.S-S. acknowledges financial support from the Spanish
I+D+i project PID2022-139555NB-I00 (TNO-JWST) funded by
MCIN/AEI/10.13039/501100011033. Based on observations made at the Sierra
Nevada Observatory (OSN), operated by the Instituto de Astrofísica de
Andalucía (IAA-CSIC), and at the Centro Astronómico Hispano-Alemán (CAHA),
operated jointly by Junta de Andalucía and the IAA-CSIC. Development of this
software would not have been possible without the invaluable and selfless
contributions of the open source community.

## References

  * Agudo et al. (2012) Agudo, I., Molina, S. N., Gómez, J. L., et al. 2012, in International Journal of Modern Physics Conference Series, Vol. 8, International Journal of Modern Physics Conference Series, 299–302, doi: 10.1142/S2010194512004746
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bokeh Development Team (2023) Bokeh Development Team. 2023, Bokeh: Python library for interactive visualization.  https://bokeh.org/
  * Bradley et al. (2023) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2023, astropy/photutils: 1.8.0, 1.8.0, Zenodo, doi: 10.5281/zenodo.7946442
  * Di Gesu et al. (2022) Di Gesu, L., Donnarumma, I., Tavecchio, F., et al. 2022, ApJ, 938, L7, doi: 10.3847/2041-8213/ac913a
  * Di Gesu et al. (2023) Di Gesu, L., Marshall, H. L., Ehlert, S. R., et al. 2023, Nature Astronomy, 7, 1245, doi: 10.1038/s41550-023-02032-7
  * Django Software Foundation (2023) Django Software Foundation. 2023, Django, 4.2.  https://www.djangoproject.com
  * Ehlert et al. (2023) Ehlert, S. R., Liodakis, I., Middei, R., et al. 2023, ApJ, 959, 61, doi: 10.3847/1538-4357/ad05c4
  * Errando et al. (2024) Errando, M., Liodakis, I., Marscher, A. P., et al. 2024, ApJ, 963, 5, doi: 10.3847/1538-4357/ad1ce4
  * Escudero Pedrosa et al. (2024) Escudero Pedrosa, J., Morcuende Parrilla, D., & Otero-Santos, J. 2024, IOP4, Zenodo, doi: 10.5281/zenodo.10222722
  * Greisen & Calabretta (2002) Greisen, E. W., & Calabretta, M. R. 2002, A&A, 395, 1061, doi: 10.1051/0004-6361:20021326
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Husillos et al. (2021) Husillos, C., Bernardos, M., Agudo, I., Bonnoli, G., & Escudero Pedrosa, J. 2021, cesarhusrod/iop3: v1.0.0, v1.0.0, Zenodo, doi: 10.5281/zenodo.10995575
  * Jester et al. (2005) Jester, S., Schneider, D. P., Richards, G. T., et al. 2005, AJ, 130, 873, doi: 10.1086/432466
  * Jordi et al. (2006) Jordi, K., Grebel, E. K., & Ammon, K. 2006, A&A, 460, 339, doi: 10.1051/0004-6361:20066082
  * Kim et al. (2024) Kim, D. E., Di Gesu, L., Liodakis, I., et al. 2024, A&A, 681, A12, doi: 10.1051/0004-6361/202347408
  * Lang et al. (2010) Lang, D., Hogg, D. W., Mierle, K., Blanton, M., & Roweis, S. 2010, AJ, 139, 1782, doi: 10.1088/0004-6256/139/5/1782
  * Liodakis et al. (2022) Liodakis, I., Marscher, A. P., Agudo, I., et al. 2022, Nature, 611, 677, doi: 10.1038/s41586-022-05338-0
  * Lupton (2005) Lupton, R. 2005, Transformations between SDSS magnitudes and other systems \- SDSS-III.  https://www.sdss3.org/dr8/algorithms/sdssUBVRITransform.php#Lupton2005
  * Marshall et al. (2023) Marshall, H. L., Liodakis, I., Marscher, A. P., et al. 2023, arXiv e-prints, arXiv:2310.11510, doi: 10.48550/arXiv.2310.11510
  * Matt Craig et al. (2022) Matt Craig, Lauren Chambers, Lauren Glattly, & Astropy Project. 2022, CCD Data Reduction Guide.  https://www.astropy.org/ccd-reduction-and-photometry-guide/
  * Middei et al. (2023a) Middei, R., Liodakis, I., Perri, M., et al. 2023a, ApJ, 942, L10, doi: 10.3847/2041-8213/aca281
  * Middei et al. (2023b) Middei, R., Perri, M., Puccetti, S., et al. 2023b, ApJ, 953, L28, doi: 10.3847/2041-8213/acec3e
  * Nilsson et al. (2007) Nilsson, K., Pasanen, M., Takalo, L. O., et al. 2007, A&A, 475, 199, doi: 10.1051/0004-6361:20077624
  * Otero-Santos et al. (2023a) Otero-Santos, J., Piirola, V., Escudero, J., et al. 2023a, The Astronomer’s Telegram, 16305, 1
  * Otero-Santos et al. (2023b) Otero-Santos, J., Piirola, V., Pacciani, L., et al. 2023b, The Astronomer’s Telegram, 16360, 1
  * Otero-Santos et al. (2024) Otero-Santos, J., Piirola, V., Escudero Pedrosa, J., et al. 2024, AJ, 167, 137, doi: 10.3847/1538-3881/ad250d
  * Patat & Romaniello (2006) Patat, F., & Romaniello, M. 2006, PASP, 118, 146, doi: 10.1086/497581
  * Peirson et al. (2023) Peirson, A. L., Negro, M., Liodakis, I., et al. 2023, ApJ, 948, L25, doi: 10.3847/2041-8213/acd242
  * Pence et al. (2010) Pence, W. D., Chiappetti, L., Page, C. G., Shaw, R. A., & Stobie, E. 2010, A&A, 524, A42, doi: 10.1051/0004-6361/201015362
  * Piirola, V. et al. (2020) Piirola, V., Berdyugin, A., Frisch, P. C., et al. 2020, A&A, 635, A46, doi: 10.1051/0004-6361/201937324
  * SQLite Developers (2023) SQLite Developers. 2023, SQLite, 3.  https://www.sqlite.org/
  * Vue Developers (2023) Vue Developers. 2023, Vue, 3.4.9.  https://vuejs.org/
  * Zapatero Osorio et al. (2005) Zapatero Osorio, M. R., Caballero, J. A., & Béjar, V. J. S. 2005, ApJ, 621, 445, doi: 10.1086/427433


lmtt

# James-Stein estimation of the first principal component

Alex Shkolnik111Department of Statistics and Applied Probability, University
of California, Santa Barbara, CA and Consortium for Data Analytics in Risk,
University of California, Berkeley, CA. Email<EMAIL_ADDRESS>
(August 31, 2021
This Version: .222I am indebted to Lisa Goldberg for conjecturing the result
derived in this work.)

###### Abstract

The Stein paradox has played an influential role in the field of high
dimensional statistics. This result warns that the sample mean, classically
regarded as the “usual estimator”, may be suboptimal in high dimensions. The
development of the James-Stein estimator, that addresses this paradox, has by
now inspired a large literature on the theme of “shrinkage” in statistics. In
this direction, we develop a James-Stein type estimator for the first
principal component of a high dimension and low sample size data set. This
estimator shrinks the usual estimator, an eigenvector of a sample covariance
matrix under a spiked covariance model, and yields superior asymptotic
guarantees. Our derivation draws a close connection to the original James-
Stein formula so that the motivation and recipe for shrinkage is intuited in a
natural way.

Keywords: Stein’s paradox, James-Stein estimator, sample eigenvectors, PCA.

## 1\. Introduction

The Stein paradox has played an influential role in the field of high
dimensional statistics. This result warns that the sample mean, classically
regarded as the “usual estimator”, may be suboptimal in high dimensions. In
particular, [stein1956] showed that the usual estimator of a location
parameter $\theta\in\mathbb{R}^{p}$ from uncorrelated Gaussian observations
becomes inadmissible when $p>2$ under a mean-squared error criterion. That is,
an estimator with a uniformly lower risk must exist. That estimator was
established by [james1961] and eponymously named.

Among the numerous perspectives that motivate the James-Stein estimator,333A
few examples include the Galtonian regression perspective promoted by
[stigler1990], the purely frequentist development of the estimator in
[gupta1991] and the geometrical explanation in [brown2012] that builds on
Stein’s original heuristic argument [stein1956, Section 1]. the empirical
Bayes perspective (see [efron1975]) is particularly elegant. Letting $\eta$
denote the sample mean computed with $n$ measurements of an unknown
$\theta\in\mathbb{R}$ and assuming an additive, normally distributed error $w$
that has a zero mean and a variance $\nu^{2}$ (e.g., $\nu=\delta/\sqrt{n}$
where each measurement has standard error $\delta$), we write

(1) $\displaystyle\eta=\theta+w\hskip 1.0pt.$

Taking a Gaussian prior on the unknown $\theta$, that is independent of $w$,
implies that

(2) $\displaystyle\textup{E}\hskip 0.5pt(\theta\hskip 1.5pt|\hskip
1.5pt\eta)=\textup{E}\hskip
0.5pt(\eta)+\bigg{(}1-\frac{\nu^{2}}{\textup{Var}(\eta)}\bigg{)}(\eta-\textup{E}\hskip
0.5pt(\eta))\hskip 1.0pt,$

the bivariate-normal conditional expectation formula. While, by definition of
conditional expectation, $\textup{E}\hskip 0.5pt(\theta\hskip 1.5pt|\hskip
1.5pt\eta)$ is the best estimator of $\theta$ in the sense of mean-squared
error, it cannot be implemented directly as the first two moments of $\eta$
are unknown.444Often, $\nu$ is assumed to be known but estimates $\hat{\nu}$
can also be used as done in [james1961]. Stein’s paradox now amounts to the
fact that “good” substitutes for $\textup{E}\hskip 0.5pt(\eta)$ and
$\textup{Var}(\eta)$ are available only in higher dimensions; precisely, when
$\theta\in\mathbb{R}^{p}$ with $p>2$.

Formula $(\ref{ce})$ extends easily to the multivariate555See formulas in
[anderson2003, Section 2.5] which can be used to design James-Stein estimators
of an unkown $\theta\in\mathbb{R}^{p}$ when $w$ has a general covariance
$\Sigma$ as in [bock1975]. case and motivates the estimator

(3) $\displaystyle\eta(\mathrm{c})=m+\mathrm{c}\hskip 1.0pt(\eta-m)$

where $m$ is an estimate (or guess) of the expected value of
$\eta\in\mathbb{R}^{p}$ and $\mathrm{c}\in(0,1)$ is a shrinkage parameter. In
words, $(\ref{jshrink})$ attempts to center the entries of $\eta$, shrinks the
resulting entries and recenters at $m$. Assuming $\nu$ is known and $p>2$,
setting

(4) $\displaystyle\mathrm{c}=1-\hskip
1.0pt\frac{\nu^{2}}{\mathrm{s}^{2}(\eta)}\Big{(}\frac{p-2}{p}\Big{)}$

where $\mathrm{s}^{2}(\eta)=\sum_{i=1}^{p}(\eta_{i}-m_{i})^{2}/p$ yields the
James-Stein estimator. Remarkably, any fixed $m\in\mathbb{R}^{p}$ (e.g.,
[stein1956] considers the origin)666 A natural choice for $m$ takes the sample
mean of $\eta$ in each entry, referred to as shrinkage toward the “grand
mean”. However, this would require $p>3$ [efron1975]).  results in an
estimator $(\ref{jshrink})$ with a strictly smaller mean-squared error than
$\eta$ [efron1975, Section 1]. While three is provably the critical
dimension777 The theme of a critical dimension is encountered frequently in
statistics and probability. [brown1971], for example, derives a close
mathematical relationship between the admissibility of the James-Stein
estimator and the transience of the Brownian motion in $\mathbb{R}^{p}$, which
also requires $p>2$., [stein1956, Section 1] heuristically argued for higher
performance in higher dimensions and for relaxing the normality.

The James-Stein estimator has inspired a rich literature on the theme of
“shrinkage” in statistics. Just a small sampling of examples includes ridge
regression [hoerl1970], the LASSO [tibshirani1996], the Ledoit-Wolf covariance
estimator [ledoit2004] and the Elastic Net [zou2005]. Excellent textbook
treatments of the ideas behind the Stein paradox and James-Stein shrinkage
include [gruber2017] and [fourdrinier2018]. In this paper, we leverage these
ideas to develop and analyze a James-Stein estimator for the first principal
component of a sample covariance matrix. The results again prove the efficacy
of James-Stein estimation, and do so for one the cornerstone methods in high-
dimensional statistics, principal component analysis [jolliffe2016].

Consider a $p\times p$ sample covariance matrix $\mathrm{S}$ that is based on
$n$ observations of some random vector $y\in\mathbb{R}^{p}$. Without loss of
generality, we write

(5) $\displaystyle\mathrm{S}=\mathscr{s}^{2}_{p}hh^{\top}+\mathrm{G}$

for $\mathrm{G}=\mathrm{S}-\mathscr{s}^{2}_{p}\hskip 1.0pthh^{\top}$ and $h$,
the sample eigenvector with the largest eigenvalue, i.e.,

(6) $\displaystyle\mathrm{S}\hskip 1.0pth=\mathscr{s}^{2}_{p}\hskip
1.0pth\hskip 16.0pt\text{and}\hskip
16.0pt\mathscr{s}^{2}_{p}=\max_{|z|=1}\langle z,\mathrm{S}z\rangle\hskip
1.0pt.$

By convention $h$ has unit length and corresponds to a direction along which
the variance of $\mathrm{S}$ (i.e., $\mathscr{s}^{2}_{p}$) is maximum (i.e.,
the first principal component). A substantial and rapidly growing literature
exists to study the ($p$ and/or $n$) asymptotic behavior of the eigenpair
$(\mathscr{s}^{2}_{p},h)$ and the remaining eigenstructure in order to
quantify either the estimation or the empirical error. See [wang2017] for
recent results and a systematic discussion of this literature. The topic of
shrinkage estimators arises naturally in this context and was raised by
[stein1986], who suggested improving the usual estimate $\mathrm{S}$ via
eigenvalue shrinkage. Indeed, there is by now a large literature on estimators
that adjust the eigenvalues of sample covariance matrices to improve their
performance with respect to some loss function [donoho2018].

In this paper, we develop and analyze a James-Stein estimator for the first
principal component of a high-dimension and low-sample (HDLS) data set.888The
HDLS framework, as introduced in [hall2005] and [ahn2007], is increasingly
relevant for data science [aoshima2018]. The recipe for the estimator begins
with $h$ and $\mathscr{s}^{2}_{p}$ in $(\ref{hands})$ and the next $(q-1)$
largest sample eigenvalues
$\mathscr{s}^{2}_{p-1},\dots,\mathscr{s}^{2}_{p-q+1}$ ($\min(n,p)>q$)
corresponding to a model with $q$ spikes.999Roughly speaking $q$ is the number
of factors (or spikes) in the data, after which a sufficiently large eigengap
(between the $q$th and the next eigenvalue) is observed (see [fan2020]).

  1. Step 1.

Set $\eta=\mathscr{s}_{p}h$, compute the sample statistics
$m(\eta)=\sum_{i=1}^{p}\eta_{i}/p$ and
$\mathrm{s}^{2}(\eta)=\sum_{i=1}^{p}(\eta_{i}-m(\eta))^{2}/p$, and define

$\displaystyle\mathrm{c}=1-\frac{\hat{\nu}^{2}}{\mathrm{s}^{2}(\eta)}\hskip
11.56346pt\text{where}\hskip
11.56346pt\hat{\nu}^{2}=\bigg{(}\frac{\textup{tr}(\mathrm{S})-(\mathscr{s}^{2}_{p}+\cdots+\mathscr{s}^{2}_{p-q+1})}{\min(n,p)-q}\bigg{)}\hskip
2.0pt/p\hskip 1.0pt.$

  2. Step 2.

Return the estimator (corrected principal component)101010With a slight abuse
of the notation, $u-x=(u_{1}-x,\dots,u_{p}-x)$ for $u\in\mathbb{R}^{p}$ and
$x\in\mathbb{R}$.

$\displaystyle h^{\hskip 1.0pt\textup{JS}}=\frac{1}{\sqrt{p}}\hskip
1.0pt\bigg{(}\frac{m(\eta)+\mathrm{c}\hskip
1.5pt(\eta-m(\eta))}{\sqrt{m^{2}(\eta)+\mathrm{c}^{2}\hskip
1.0pt\mathrm{s}^{2}(\eta)}}\bigg{)}\hskip 1.0pt.$

The vector $h^{\hskip 1.0pt\textup{JS}}$ is the James-Stein estimator of the
first principal component of the data. The numerator contains the shrinkage
formula $(\ref{jshrink})$ while the divisor normalizes the shrunk vector to a
unit length (by convention). The relationship to the shrinkage parameter
$\mathrm{c}$ in $(\ref{jsc})$ is evident by treating $p$ as large. The
estimate $\hat{\nu}^{2}$ corresponds to the bulk of the eigenvalue spectrum,
and may be viewed the “noise” in the context of a signal-to-noise ratio that
plays a prominent role of the results in Sections 2 & 3.

It is reasonable to suspect that a James-Stein type shrinkage of the principal
component $h$, a high dimensional vector, could improve either the convergence
rate or accuracy of the limit in some appropriate asymptotic regime. However,
the standard orthonormal transformation and the eigengap partition of the
sample eigenvectors, that is typically leveraged by their asymptotic analyses
(e.g., [paul2007], [shen2016] and [wang2017]), can obscure the systematic
nature of the sample bias. As sensibly pointed out by [wang2017] in reference
to the partition of the sample eigenvector, the “two parts intertwine in such
a way that correction for the biases of estimating eigenvectors is almost
impossible.” However, in the original (untransformed) coordinate system and
the HDLS asymptotic regime, the bias can in fact be identified, characterized
and (partially) corrected. This program was carried out in [goldberg2021], who
adopt a factor model in an HDLS regime and utilize a portfolio theory
application to motivate their analysis.

The main results of this paper rederive the adjustment of [goldberg2021] but
within a James-Stein type framework. In particular, we establish identity
$(\ref{jsprob})$ in which $\eta=\mathscr{s}_{p}h$ and $\theta$ is related to
the associated population eigenvector. From here, the James-Stein shrinkage
acts on the perturbation $w$ so that the estimator $h^{\hskip
1.0pt\textup{JS}}$ outperforms $h$ on the mean-squared error and angle metrics
as $p\uparrow\infty$. The theoretical guarantees provided here are new and
their proofs rely on a different set of mathematical tools than
[goldberg2021]. In particular, the new approach leverages Weyl’s inequality
and the Davis-Kahan theorem from matrix perturbation theory to give simpler
proofs and potentially expand the scope of applicability of the resulting
estimator. The HDLS regime, in which the number of variables $p$ grows to
infinity, and the number of observations $n$ to be fixed, plays a crucial role
in the analysis.

The paper is organized as follows. Section 2 defines the spiked covariance
model underlying our results. Section 3 develops the James-Stein estimator
$h^{\hskip 1.0pt\textup{JS}}$ and Section 4 proves the theoretical guarantees
for this estimator. Appendix A contains proofs of the auxiliary results. The
following notation is used throughout. Let $\langle u,v\rangle$ denote the
standard inner product of $u,v\in\mathbb{R}^{d}$ so that $|u|=\sqrt{\langle
u,u\rangle}$ and $m(u)=\langle u,\mathrm{e}\rangle/d$ where
$\mathrm{e}=(1,\dots,1)^{\top}$ are the length and mean. Set
$\mathrm{s}^{2}(u)=|u-m(u)|^{2}/d$ and $\textup{cov}(u,v)=\langle
u-m(u),v-m(v)\rangle/d$ (see the notation of footnote 10). We use a subscript
$1\leq p\leq\infty$ to highlight the dependence on $p$ of various quantities,
e.g., $m(\eta)=m_{p}(\eta)$ for $\eta\in\mathbb{R}^{p}$ and $m_{\infty}(\eta)$
is the limit $\lim_{p\uparrow\infty}m_{p}(\eta)$ when it exists.

## 2\. A scarcely sampled spiked model

We use a spiked covariance model, borrowed from the HDLS literature. We also
restrict ourselves to a single unbounded spike in the “boundary case” (see
[jung2012]) wherein the largest eigenvalue of the covariance matrix grows
linearly in $p$. In particular, consider a mean-zero (w.l.o.g.) random vector
$y\in\mathbb{R}^{p}$ with a $p\times p$ covariance matrix
$\Sigma=\textup{Var}(y)$ and let,

(7) $\displaystyle\Sigma=\Gamma+\beta\beta^{\top}$

for a symmetric, positive-semidefinite $p\times p$ matrix $\Gamma$ and vector
$\beta\in\mathbb{R}^{p}$. The following affirms that $b=\beta/|\beta|$ is an
eigenvector of $\Sigma$ with eigenvalue $\langle\beta,\beta\rangle$.

######  Assumption 2.1 (w.l.o.g.).

$\Gamma b=0$ and $m_{p}(b)\geq 0$ for any $p$.

To state our additional assumptions on the model, we project the data vector
$y\in\mathbb{R}^{p}$ onto the eigenvectors of $\Sigma$. More precisely, define

(8)
$\displaystyle\psi_{p}=\langle\beta,y\rangle/\langle\beta,\beta\rangle\hskip
1.0pt.$

It is immediate that $\textup{E}\hskip 0.5pt(\psi_{p})=0$ and
$\textup{Var}(\psi_{p})=1$. For every eigenvalue $\alpha_{i}$ of $\Gamma$, let

(9) $\displaystyle\phi_{i}=\langle\gamma,y\rangle\hskip
16.0pt\text{where}\hskip 16.0pt\gamma\in\mathbb{R}^{p}\hskip 2.0pt\hskip
0.5pt:\hskip 0.5pt\hskip 2.0pt\Gamma\hskip 1.0pt\gamma=\alpha_{i}\hskip
1.0pt\gamma\hskip 1.0pt.$

Note, $\textup{E}\hskip 0.5pt(\phi_{i})=0$ and
$\textup{Var}(\phi_{i})=\alpha_{i}$. As the dimension $p$ grows we obtain a
sequence $\\{\phi_{i}\\}_{i\geq 1}$. As a technical remark,
$\phi_{i}=\phi_{i,p}$ and $\alpha_{i}=\alpha_{i,p}$ depend on $p$. We consider
a sequence of models $(\ref{popcov})$ constructed from sequences
$\\{\beta_{i}\\}_{i\geq 1}$ and $\\{\Gamma_{p}\\}$.

######  Assumption 2.2.

For constants $\mu\in\mathbb{R}$ and $\sigma,\delta\in(0,\infty)$ as
$p\uparrow\infty$ we have:

  1. (i)

$m_{\infty}(\beta)=\mu$ and $\mathrm{s}^{2}_{\infty}(\beta)=\sigma^{2}$.

  2. (ii)

$\psi_{\infty}=\lim_{p\uparrow\infty}\psi_{p}$ exists as a $\mathbb{R}$-valued
random variable almost surely.

  3. (iii)

$m_{\infty}(\phi)=m_{\infty}(\varphi)=0$ almost surely for
$\\{\varphi_{i}\\}_{i\geq 1}$ with
$\varphi_{i}=\varphi_{i,p}=\phi_{i}pm(\gamma^{i})$.

  4. (iv)

$\mathrm{s}^{2}_{\infty}(\phi)=m_{\infty}(\alpha)=\delta^{2}$ almost surely.

Condition (i) imposes regularity on the sequence $\\{\beta_{i}\\}_{i\geq 1}$
and implies that the largest eigenvalue of $\Sigma$ (i.e.,
$\langle\beta,\beta\rangle$) grows linearly with $p$. The random variable
$\psi_{\infty}$ in (ii) is closely related to a principal component score (in
the limit $p\uparrow\infty$), and it captures the randomness along the first
(population) principal component. Conditions (iii) and (iv) are related to
certain requirements on a measure of sphericity of the model (e.g.,
$\textup{tr}(\Sigma)^{2}/(\textup{tr}(\Sigma^{2})p)$) and summarize the
conclusions of [jung2009, Theorem 1]. In particular, (iii) may be viewed as
laws of large numbers for $\\{\phi_{i}\\}_{i\geq 1}$ and
$\\{\varphi_{i}\\}_{i\geq 1}$ and suggests $pm(\gamma^{i})\asymp 1$ whereas
(i) implies the first principal component has $\sqrt{p}m(b)\asymp 1$. This
highlights that the spike eigenvector $b$ differs from those of $\Gamma$ both
in terms of the magnitude of its eigenvalue as well as the structure of the
vector itself. According to (iv), the average eigenvalue of $\Gamma$ (unlike
the spike eigenvalue $\langle\beta,\beta\rangle$) is bounded in $p$, and per
(iii), the eigenvectors $\gamma^{i}$ do not have entries biased towards a
nonzero mean $\sqrt{p}m(\gamma^{i})$, unlike $b$ (an exception is $\mu=0$,
which is the case when the James-Stein estimator will turn out to be
ineffective; see Remark 2.8).

The following assumption is a standard one in statistical data analysis but
may be (and often is) relaxed in the HDLS setup (e.g., [jung2012]).

######  Assumption 2.3.

There are a fixed $n\geq 2$ i.i.d. observations of $y\in\mathbb{R}^{p}$.

Our forthcoming results hold even when only $2$ observations are available,
hence, a scarcely sampled model. Let $\mathrm{Y}$ be the $p\times n$ data
matrix with the $k$th column containing the $k$th observation of $y$, and
define the sample covariance matrix $\mathrm{S}$ by

(10) $\displaystyle\mathrm{S}=\mathrm{Y}\mathrm{Y}^{\top}/n\hskip 1.0pt.$

We let $h\in\mathbb{R}^{p}$ denote the eigenvector of $\mathrm{S}$ with the
largest eigenvalue $\mathscr{s}^{2}_{p}$ (see $(\ref{hands})$). It is unique
only up to sign (and $|h|=1$), motivating the following (c.f. Assumption 2.1).

######  Assumption 2.4 (w.l.o.g.).

$m_{p}(h)\geq 0$ for any $p$.

We write $\mathrm{S}=\mathrm{G}+\eta\eta^{\top}$ in analogy to
$(\ref{popcov})$ (taking $\mathrm{G}=\mathrm{S}-\eta\eta^{\top}$) and set

(11) $\displaystyle\eta=\mathscr{s}_{p}h\hskip 1.0pt.$

Next, define the following measure of finite-sample distortion,

(12) $\displaystyle\chi_{n}^{2}=\langle\mathcal{X},\mathcal{X}\rangle/n\hskip
16.0pt\text{and}\hskip
16.0pt\mathcal{X}=\lim_{p\uparrow\infty}\frac{\mathrm{Y}^{\top}\beta}{\langle\beta,\beta\rangle}\hskip
1.0pt.$

The latter limit exists by Assumption 2.2 while Assumption 2.3 implies that
$\mathcal{X}\in\mathbb{R}^{n}$ has i.i.d. entries (distributed as
$\psi_{\infty}$). Consequently, we have $\chi^{2}_{n}\to 1$ as
$n\uparrow\infty$.

We can measure the error in any estimator $\eta$ of $\theta$ by the mean-
squared error, as would be consistent with the James-Stein framework.

(13) $\displaystyle\textup{MSE}_{p}(\eta\hskip 1.5pt|\hskip
1.5pt\theta)=\langle\eta-\theta,\eta-\theta\rangle/p$

######  Proposition 2.5.

Let $\theta=\chi_{n}\beta$ and suppose Assumptions 2.1–2.4 hold. Then,

(14) $\displaystyle\textup{MSE}_{\infty}(\eta\hskip 1.5pt|\hskip
1.5pt\theta)=\frac{\hskip 1.0pt\delta^{2}}{n\hskip 1.0pt}\hskip 1.0pt.$

######  Remark 2.6.

If $\theta=\beta$ then the right side would be multiplied by the factor,
$1+\frac{\textup{SNR}^{2}}{\mathrm{r}^{2}_{\infty}}\hskip
1.0pt\Big{(}\frac{\chi_{n}-1}{\chi_{n}}\Big{)}^{2}$ where we define a signal-
to-noise ratio SNR and signal-incoherence $\mathrm{r}_{\infty}$ as

(15) $\displaystyle\textup{SNR}=\Big{(}\frac{\sigma}{\delta}\Big{)}\hskip
1.0pt\chi_{n}\sqrt{n}\hskip 11.56346pt\text{and}\hskip
11.56346pt\mathrm{r}_{\infty}=\frac{1}{\sqrt{1+(\mu/\sigma)^{2}}}\hskip
1.0pt.$

For SNR, we regard $\sigma\chi_{n}$ as a distorted signal, and
$\delta/\sqrt{n}$ as noise that vanishes as $n$ grows (also,
$\sigma\chi_{n}\to\sigma$). The signal-incoherence $\mathrm{r}_{\infty}$ is
the limit of
$\mathrm{r}_{p}=\mathrm{r}_{p}(\beta)=\mathrm{s}_{p}(\beta)/|\beta|$ (per
Assumption 2.2) determined by the signal-to-noise ratio $\mu/\sigma$ of the
vector $\beta$. A large value of $\mathrm{r}_{\infty}$ corresponds to more
variation, or “incoherence” in the entries $\\{\beta_{i}\\}_{i\geq 1}$

A more standard way to evaluate the goodness of a sample eigenvector $h$ is
via its angle away from its population counterpart $b=\beta/|\beta|$. To this
end, let

(16) $\displaystyle\textup{SPH}_{p}(h\hskip 1.5pt|\hskip
1.5ptb)=\textup{SPH}_{p}(\eta\hskip 1.5pt|\hskip
1.5pt\theta)=\sin^{2}\Big{(}\arccos\frac{\langle\eta,\theta\rangle}{|\eta||\theta|}\Big{)}\hskip
1.0pt.$

######  Proposition 2.7.

Let $\theta=\chi_{n}\beta$ and suppose Assumptions 2.1–2.4 hold. Then,

(17) $\displaystyle\textup{SPH}_{\infty}(\eta\hskip 1.5pt|\hskip
1.5pt\theta)=\frac{\mathrm{r}^{2}_{\infty}}{\mathrm{r}^{2}_{\infty}+\textup{SNR}^{2}}\hskip
2.0pt.$

######  Remark 2.8.

Clearly, $(\ref{spherr})$ also holds with $\theta=\beta$.

## 3\. James-Stein estimation of sample eigenvectors

Having established that the sample eigenvector corresponding to the spike
(i.e., the largest eigenvalue) carries finite sample error, it is natural to
ask whether James-Stein shrinkage can improve this “usual” estimator. The key
to this question is the (to be established) identity

(18) $\displaystyle\eta=\theta+w\hskip 16.0pt\text{and}\hskip
16.0pt\theta=\chi_{n}\beta$

for a random vector $w\in\mathbb{R}^{p}$ specified in $(\ref{dev})$ of Section
4.1. As in definition $(\ref{etah})$, we have $\eta=\mathscr{s}_{p}h$ where
$h$ is the sample eigenvector with the largest eigenvalue,
$\mathscr{s}^{2}_{p}$. The perturbation $w$ of $\theta$ turns out to be such
that the shrinkage of $\eta$ is effective. We remark that the recipe of
Section 1 extends our derivation of the James-Stein estimator below to the
case of multiple (there, $q$) spikes in a natural way. This extension is
effective because the eigenvectors corresponding to the spikes are mutually
orthogonal, but it is suboptimal. An optimal estimator in the multi-spiked
setup is left for future work.

### 3.1. The JS estimator

Equation $(\ref{jsform})$ establishes a relationship between the sample and
population eigenvectors that suggests a James-Stein estimator may be derived.
An informal derivation proceeds as follows. Consider the shrinkage parameter

(19) $\displaystyle\mathrm{c}=1-\frac{\hat{\nu}^{2}}{\mathrm{s}^{2}(\eta)}$

based on $(\ref{jsc})$ with $\hat{\nu}$, an estimate of the “noise”. It is
reasonable to assign the latter to be the average of the non-spiked, non-zero
eigenvalues of $\mathrm{S}$. That is,

(20) $\displaystyle\hskip
92.50551pt\hat{\nu}^{2}=\bigg{(}\frac{\textup{tr}(\mathrm{S})-\mathscr{s}^{2}_{p}}{n-1}\bigg{)}\hskip
2.0pt/p\hskip 46.25275pt(p\geq n),$

where the scaling by $p$ turns out to be necessary due to the counterintuitive
behavior of the HDLS asymptotics. When $p<n$ the divisor $n-1$ must be
replaced by $p-1$.

This paves the way for the James-Stein sample eigenvector estimate,

$\displaystyle\eta^{\hskip 1.0pt\textup{JS}}=m(\eta)+\mathrm{c}\hskip
2.0pt(\eta-m(\eta))$

of the unnormalized eigenvector and by convention, we take unit length
version,

(21) $\displaystyle h^{\hskip 1.0pt\textup{JS}}=\frac{\eta^{\hskip
1.0pt\textup{JS}}}{|\eta^{\hskip 1.0pt\textup{JS}}|}=\frac{1}{\sqrt{p}}\hskip
1.0pt\bigg{(}\frac{m(h)+\mathrm{c}\hskip
1.5pt(h-m(h))}{\sqrt{m^{2}(h)+\mathrm{c}^{2}\hskip
1.0pt\mathrm{s}^{2}(h)}}\bigg{)}$

as the James-Stein estimator of the population eigenvector $b=\beta/|\beta|$.

The following James-Stein type theorems characterize the improvement due to
shrinkage in the original mean-squared sense as well as in the angle metrics.

######  Theorem 3.1.

Suppose Assumptions 2.1–2.4 hold. Then, almost surely,

$\displaystyle\textup{MSE}_{\infty}(\eta^{\hskip 1.0pt\textup{JS}}\hskip
1.5pt|\hskip 1.5pt\theta)$ $\displaystyle=\mathrm{c}_{\infty}\hskip
1.0pt\textup{MSE}_{\infty}(\eta\hskip 1.5pt|\hskip 1.5pt\theta)$

where $\mathrm{c}_{\infty}\in(0,1)$ is the limit of
$\mathrm{c}_{p}=\mathrm{c}$ in $(\ref{jsceig})$ with SNR defined in
$(\ref{snr})$ and

(22)
$\displaystyle\mathrm{c}_{\infty}=\frac{\textup{SNR}^{2}}{1+\textup{SNR}^{2}}\hskip
1.0pt.$

######  Theorem 3.2.

Suppose Assumptions 2.1–2.4 hold. Then, almost surely,

$\displaystyle\textup{SPH}_{\infty}(\eta^{\hskip 1.0pt\textup{JS}}\hskip
1.5pt|\hskip 1.5pt\theta)=\textup{SPH}_{\infty}(h^{\hskip
1.0pt\textup{JS}}\hskip 1.5pt|\hskip 1.5ptb)$
$\displaystyle=\mathrm{d}_{\infty}\hskip 1.0pt\textup{SPH}_{\infty}(h\hskip
1.5pt|\hskip 1.5ptb)$

where $\mathrm{d}_{\infty}\in[\mathrm{c}_{\infty},1]$ where
$\mathrm{c}_{\infty}$ is in $(\ref{clim})$ and with SNR and
$\mathrm{r}_{\infty}$ in $(\ref{snr})$, we have

(23)
$\displaystyle\mathrm{d}_{\infty}=\mathrm{c}_{\infty}+\frac{\mathrm{r}^{2}_{\infty}}{1+\textup{SNR}^{2}}\hskip
1.0pt.$

Related results may be found in [goldberg2020] and [goldberg2021] but the
metrics there are motivated by solutions of certain quadratic programs that
are useful in finance and portfolio theory.

######  Remark 3.3.

Note, $\mathrm{d}_{\infty}=1$ (i.e., no improvement in angle) if and only if
$\mu=0$.

### 3.2. The geometry of Stein’s paradox

We shed insight into the James-Stein estimator in $(\ref{hJS})$ by deriving
general conditions under which Theorems 3.1 and 3.2 hold. Our analysis adopts
the pure frequentist perspective in [gupta1991] and supplements it by
illustrating the Euclidean and the spherical geometry of the estimator. The
two geometries reflect the definitions of the error (MSE and SPH) in the two
theorems.

######  Lemma 3.4.

Let $\eta=\theta+w$ for $\theta,w\in\mathbb{R}^{p}$. Then, the solutions of
the optimizations $\min_{c\in\mathbb{R}}\textup{MSE}\hskip 1.0pt(\eta(c)\hskip
1.5pt|\hskip 1.5pt\theta)$ and $\min_{c\in\mathbb{R}}\textup{SPH}\hskip
1.0pt(\eta(c)\hskip 1.5pt|\hskip 1.5pt\theta)$ (see $(\ref{mse})$ and
$(\ref{sph})$) are given by

(24)
$\displaystyle\mathrm{c}^{\textup{MSE}}=\frac{\textup{cov}(\theta,\eta)}{\mathrm{s}^{2}(\eta)}\hskip
16.0pt\text{and}\hskip
16.0pt\mathrm{c}^{\textup{SPH}}=\frac{m(\eta)}{m(\theta)}\hskip
2.0pt\mathrm{c}^{\textup{MSE}}\hskip 1.0pt.$

The next assumptions may be viewed as laws of large numbers in the random
setting or regularity conditions in a deterministic one. They concern the
sequences $\\{\theta_{i}\\}_{i=1}^{\infty}$ and $\\{w_{i}\\}_{i=1}^{\infty}$,
and allow for dependence on $p$ (i.e., $\theta_{i}=\theta_{i}^{(p)}$ and
$w_{i}=w_{i}^{(p)}$).

######  Assumption 3.5.

For constants $m\in\mathbb{R}$ and $\nu,\xi\in(0,\infty)$ as
$p\uparrow\infty$, we have:

  1. (i)

$m_{\infty}(\theta)=m$ and $\mathrm{s}^{2}_{\infty}(\theta)=\xi^{2}$,

  2. (ii)

$m_{\infty}(w)=0$ and $\mathrm{s}^{2}_{\infty}(w)=\nu^{2}$,

  3. (iii)

$\textup{cov}_{\infty}(\theta,w)=0\hskip 1.0pt$,

  4. (iv)

there exists an estimator $\hat{\nu}=\hat{\nu}_{p}$ for each $p$ with
$\hat{\nu}_{\infty}=\nu$.

The following identities follow by direct calculation.

######  Lemma 3.6.

Suppose $\\{\theta_{i}\\}$ and $\\{w_{i}\\}$ satisfy Assumption 3.5 and
$\eta_{i}=\theta_{i}+w_{i}$. Then (almost surely), $m_{\infty}(\eta)=m$,
$\textup{cov}_{\infty}(\eta,\theta)=\xi^{2}$ and
$\mathrm{s}^{2}_{\infty}(\eta)=\xi^{2}+\nu^{2}$.

We define the signal-to-noise ratio SNR and the signal-incoherence
$\mathrm{r}_{\infty}$ as

(25) $\displaystyle\textup{SNR}=\frac{\xi}{\nu}\hskip
11.56346pt\text{and}\hskip
11.56346pt\mathrm{r}_{\infty}=\frac{1}{\sqrt{1+(m/\xi)^{2}}}\hskip 1.0pt,$

which are compatible with $(\ref{snr})$ upon taking $\theta=\chi_{n}\beta$ and
$\nu=\delta/\sqrt{n}$. The following result establishes the conclusions of
Theorems 3.1 and 3.2 in our abstract setting.

######  Proposition 3.7.

Let $\eta=\theta+w$ where $\theta,w\in\mathbb{R}^{p}$ and an estimator
$\hat{\nu}$ satisfy Assumption 3.5. Then, for $\mathrm{c}_{\infty}$ and
$\mathrm{d}_{\infty}$ defined in $(\ref{clim})$ and $(\ref{slim})$ but with
SNR and $\mathrm{r}_{\infty}$ in $(\ref{asnr})$ the estimate
$\eta(\mathrm{c})=\eta+\mathrm{c}\hskip 1.0pt(\eta-m(\eta))$ with parameter
$\mathrm{c}=1-\frac{\hat{\nu}^{2}}{\mathrm{s}^{2}(\eta)}$ satisfies

$\displaystyle\textup{MSE}_{\infty}(\eta(\mathrm{c})\hskip 1.5pt|\hskip
1.5pt\theta)=\mathrm{c}_{\infty}\hskip 1.0pt\textup{MSE}_{\infty}(\eta\hskip
1.5pt|\hskip 1.5pt\theta)\hskip 5.78172pt\text{and}\hskip
5.78172pt\textup{SPH}_{\infty}(\eta(\mathrm{c})\hskip 1.5pt|\hskip
1.5pt\theta)=\mathrm{d}_{\infty}\hskip 1.0pt\textup{SPH}_{\infty}(\eta\hskip
1.5pt|\hskip 1.5pt\theta)\hskip 1.0pt.$

Moreover, the optimal parameters $\mathrm{c}^{\textup{MSE}}$ and
$\mathrm{c}^{\textup{SPH}}$ in $(\ref{cex})$ converge as $p\uparrow\infty$ to
$\mathrm{c}_{\infty}$.

Figure 1 illustrates the geometry of the estimator $\eta(\mathrm{c})$ of the
vector $\theta$.

Figure 1: Illustration of the estimator $\eta(\mathrm{c})$ in low _(top left)_
, high _(top right)_ , and limiting _(bottom left)_ dimensions, relative to
the shrinkage target $m\in\mathbb{R}^{p}$, the vector with all entries equal
to $m(\eta)$. The open circle marks the estimator with the optimal shrinkage
parameter $\mathrm{c}^{\textup{MSE}}$. In Euclidean geometry, the estimator
$\eta(\mathrm{c})$ is located anywhere on the ray originating at the
observation $\eta$ and passing through $m$, because $\mathrm{c}\leq 1$ while
$\eta(1)=\eta$ and $\eta(0)=m$. The spherical geometry _(bottom right)_
presents the limiting ($p=\infty$) analog of the illustration in the bottom
left. There, $b=\theta/|\theta|$, $h=\eta/|\eta|$ and $z=m/|m|$ and the
contour describes all vectors with mean entry $m(b)$, which highlights the
difference between the two geometries.

## 4\. Proofs the main results

We proceed by the following three main steps.

  1. (1)

We establish the key identity $\eta=\theta+w$ with $\theta=\chi_{n}\beta$ per
$(\ref{jsform})$ in Section 4.1.

  2. (2)

We derive the convergence properties in the HDLS regime of the eigenvalues and
eigenvectors of $\mathrm{S}$ under our spiked covariance model setting in
Section 4.2.

  3. (3)

We verify that $\theta,w$ in (1) and $\hat{\nu}$ in $(\ref{noisest})$ satisfy
the conditions of Assumption 3.5, which leads to the guarantees for James-
Stein shrinkage in Proposition 3.7.

Theorems 3.1 and 3.2 and then corollaries of Proposition 3.7, which is proved
in Appendix A. We will make use of two classic results in matrix perturbation
theory.

######  Theorem (Weyl).

Let $\mathrm{A}$ and $(\mathrm{A}+\Delta)$ be (real) symmetric $n\times n$
matrices with eigenvalues $\alpha_{1}\geq\cdots\geq\alpha_{n}$ and
$\zeta_{1}\geq\cdots\geq\zeta_{n}$ respectively. Then,

$\displaystyle\max_{1\leq j\leq n}|\alpha_{j}-\zeta_{j}|\leq|\Delta|\hskip
1.0pt.$

For a proof see [horn2012] (also [weyl1912]).

######  Theorem (Davis-Kahan).

Let $\mathrm{A}$ and $(\mathrm{A}+\Delta)$ be (real) symmetric $n\times n$
matrices with $\mathrm{A}\hskip 1.0pta^{j}=\alpha_{j}a^{j}$ and
$(\mathrm{A}+\Delta)\hskip 1.0ptb^{j}=\beta_{j}b^{j}$ for eigenvectors
$a^{j},b^{j}\in\mathbb{R}^{n}$ and eigenvalues
$\alpha_{j},\beta_{j}\in\mathbb{R}$. Suppose
$\alpha_{1}\geq\cdots\geq\alpha_{n}$ and $\beta_{1}\geq\cdots\geq\beta_{n}$
with the convention $\alpha_{0}=\infty=-\alpha_{n+1}$ and assume
$\gamma_{j}=\min\\{\alpha_{j-1}-\alpha_{j},\alpha_{j}-\alpha_{j+1}\\}>0$.
Then,

$\displaystyle\hskip 23.12692pt|\hskip 1.0pta^{j}-b^{j}|\leq\frac{3}{\hskip
2.0pt\gamma_{j}}\hskip 1.0pt|\Delta|\hskip 23.12692pt\text{provided (w.l.o.g.)
}\langle a^{j},b^{j}\rangle\geq 0.$

This result is proved in [yu2015, Corollary 1].

### 4.1. Establishing the key identity

A key tool for random matrix theory in the HDLS regime is the dual sample
covariance matrix. This is $n\times n$ matrix ($n\geq 2$ is fixed),

(26) $\displaystyle\mathrm{L}=\mathrm{Y}^{\top}\mathrm{Y}/p\hskip 1.0pt.$

The next result is well known and relates the spectra of
$\mathrm{S}=\mathrm{Y}\mathrm{Y}^{\top}/n$ and $\mathrm{L}$.

######  Lemma 4.1.

Let $\mathrm{L}u=\ell^{2}u$ where $\ell^{2}\in(0,\infty)$ and
$u\in\mathbb{R}^{n}$. Then, $\mathrm{S}v=\mathscr{s}^{2}v$ where
$v=\mathrm{Y}u/(\sqrt{p}\ell)$ and $\mathscr{s}^{2}=\ell^{2}p/n$. Conversely,
let $\mathrm{S}v=\mathscr{s}^{2}v$ where $\mathscr{s}^{2}\in(0,\infty)$ and
$v\in\mathbb{R}^{p}$. Then, $\mathrm{L}u=\ell^{2}u$ where
$u=\mathrm{Y}^{\top}v/(\sqrt{n}\mathscr{s})$ and
$\ell^{2}=\mathscr{s}^{2}n/p$.

Proof. Multiplying the identity $\mathrm{L}u=\ell^{2}u$ by $\mathrm{Y}$ from
both sides, we obtain

$\displaystyle\mathrm{Y}\mathrm{L}u=\ell^{2}\mathrm{Y}u\hskip
11.56346pt\Rightarrow\hskip
11.56346pt\mathrm{S}\mathrm{Y}u=\Big{(}\frac{\ell^{2}p}{n}\Big{)}\mathrm{Y}u\hskip
1.0pt.$

Note, $(\mathrm{Y}u)^{\top}(\mathrm{Y}u)=(u^{\top}\mathrm{L}u)p=\ell^{2}p$, so
$v=(\mathrm{Y}u)/(\sqrt{p}\ell)$ has unit length. Dividing by $\sqrt{p}\ell$
yields $\mathrm{S}v=\mathscr{s}^{2}v$ as required. The converse has an
identical argument. ∎

The spike model $\Sigma=\Gamma+\beta\beta^{\top}$ has a full basis of
eigenvectors given by $b=\beta/|\beta|$ and $\\{\gamma^{i}\\}_{i=1}^{p-1}$,
the latter corresponding to the nonzero eigenvalues $\\{\alpha\\}_{i=1}^{p-1}$
of $\Gamma$. Thus,

$\displaystyle y$ $\displaystyle=\langle b,y\rangle
b+{\textstyle\sum\nolimits}_{i=1}^{p}\langle\gamma^{i},y\rangle\hskip
1.0pt\gamma^{i}=\beta\psi+\epsilon$

where
$\epsilon={\textstyle\sum\nolimits}_{i=1}^{p}\phi_{i}\gamma^{i}\in\mathbb{R}^{p}$
and $\psi=\psi_{p}=\langle\beta,y\rangle/\langle\beta,\beta\rangle$ as in
$(\ref{xp})$. Consequently, letting $\mathrm{Y}$ denote the $p\times n$ matrix
of i.i.d. observations of $y\in\mathbb{R}^{p}$, we have

$\displaystyle\mathrm{Y}=\beta\mathrm{X}^{\top}+\mathcal{E}\hskip 1.0pt.$

where
$\mathrm{X}=\mathrm{Y}^{\top}\beta/\langle\beta,\beta\rangle\in\mathbb{R}^{n}$
consists of i.i.d. observations of $\psi$ and $\mathcal{E}$ is a $p\times n$
matrix with i.i.d. columns consisting of the observations of $\epsilon$ as
defined above.

By orthogonality, we obtain that
$\mathrm{L}=\mathrm{Y}^{\top}\mathrm{Y}/p=(\langle\beta,\beta\rangle/p)\hskip
2.0pt\mathrm{X}\mathrm{X}^{\top}+\mathcal{E}\mathcal{E}^{\top}/p$. Let $x_{p}$
be the eigenvector of $\mathrm{L}$ with the largest eigenvalue,
$\ell_{p}^{2}$, and $x_{\infty}=\mathcal{X}/\chi_{n}\in\mathbb{R}^{n}$ (the
unit length normalization of $\mathcal{X}$), where $\mathcal{X}$ and
$\chi_{n}$ are defined in $(\ref{dst})$. By Lemma 4.1,

$\displaystyle h$
$\displaystyle=\frac{\mathrm{Y}x_{p}}{\sqrt{p}\ell_{p}}=\frac{(\beta\mathrm{X}^{\top}+\mathcal{E})x}{\sqrt{p}\ell_{p}}=\beta\hskip
1.0pt\frac{\chi_{n}\langle
x_{p},x_{\infty}\rangle}{\sqrt{p/n}\ell_{p}}+\frac{\mathcal{E}x_{p}}{\sqrt{p}\ell_{p}}$
$\displaystyle=\beta\Big{(}\frac{\chi_{n}}{\mathscr{s}_{p}}\Big{)}+\beta\Big{(}\frac{\chi_{n}}{\mathscr{s}_{p}}\Big{)}\langle
x_{p},x_{p}-x_{\infty}\rangle+\frac{\mathcal{E}x_{p}}{\sqrt{n}\mathscr{s}_{p}}$

We deduce that $\eta=\mathscr{s}_{p}h=\chi_{n}\beta+w$ as required by
$(\ref{jsform})$ where

(27) $\displaystyle w=\chi_{n}\beta\hskip 1.0pt\langle
x_{p},x_{p}-x_{\infty}\rangle+\frac{\mathcal{E}x_{p}}{\sqrt{n}}\hskip 1.0pt.$

### 4.2. Convergence of the eigen- values/vectors

It is not difficult to establish that limit of $\mathrm{L}=\mathrm{L}^{(p)}$
as $p\uparrow\infty$ (in any norm on $\mathbb{R}^{n\times n}$) takes the
following form

$\displaystyle\mathrm{L}^{(\infty)}=(\sigma^{2}+\mu^{2})\hskip
1.0pt(n\chi^{2}_{n})\hskip
2.0ptx_{\infty}x_{\infty}^{\top}+\delta^{2}\mathrm{I}\hskip 1.0pt.$

The first term is the limit of $(\langle\beta,\beta\rangle/p)\hskip
1.0pt\mathrm{X}\mathrm{X}^{\top}$ under Assumption 2.2 (note that
$\langle\beta,\beta\rangle=\mathrm{s}^{2}(\beta)+m^{2}(\beta)$) and the
definitions of $\mathrm{X}$ and $x_{\infty}$ above. By Assumption 2.3, the
columns of $\mathcal{E}$ are i.i.d. copies of $\epsilon$ with
$\textup{E}\hskip 0.5pt(\epsilon)=0_{p}$ by definition, the strong law of
large numbers, confirms the that the off-diagonal entries of the second term
are zero. That $\delta^{2}$ determines all the diagonal entries is again a
consequence of Assumption 2.2. This entails proving that
$\langle\epsilon,\epsilon\rangle/p$ converges almost surely to $\delta^{2}$ as
is done in Section 4.3 item (ii).

The matrix $\mathrm{L}^{(\infty)}$ has an easily described spectrum. Its
largest eigenvalue is given by $\ell_{\infty}^{2}=(\sigma^{2}+\mu^{2})\hskip
1.0pt(n\chi^{2}_{n})+\delta^{2}$ and has the eigenvector $x_{\infty}$. All
remaining eigenvalues equal $\delta^{2}$. Since $\mathrm{L}^{(p)}$ converges
(in any norm) to $\mathrm{L}^{(\infty)}$, the largest eigenvalue
$\ell^{2}_{p}$ converges to $\ell^{2}_{\infty}$ almost surely. All the
remaining eigenvalues converge to $\delta^{2}$.

By Weyl’s inequality, setting $\mathrm{A}=\Omega$ and
$\mathrm{B}=(\mathrm{L}-\Omega)$ so that $\mathrm{A}+\mathrm{B}=\mathrm{L}$,
we have

(28) $\displaystyle|\ell^{2}_{p-i+1}-((\sigma^{2}+\mu^{2})\hskip
1.0pt(n\chi^{2}_{n})1_{\\{i=1\\}}+\delta^{2})|$
$\displaystyle\leq|\mathrm{L}^{(p)}-\mathrm{L}^{(\infty)}|\hskip 16.0pt1\leq
i\leq n\hskip 1.0pt.$

We immediately deduce (by Lemma 4.1) the following result, variants of which
appear in [jung2012], [shen2016] and [goldberg2021]. Let
$\mathscr{s}^{2}_{p-i+1}$ denote the $i$th largest eigenvalue of $\mathrm{S}$
(for $i>\min(p,n)$ all eigenvalues are zero).

######  Proposition 4.2.

Fix $n\geq 2$ and suppose Assumptions 2.1, 2.2 and 2.3 hold. Then,
$\lim_{p\uparrow\infty}\mathscr{s}^{2}_{p-i+1}/p=1_{\\{i=1\\}}\chi_{n}^{2}(\sigma^{2}+\mu^{2})+\delta^{2}/n$
almost surely for fixed $1\leq i\leq n$.

Next, by the Davis-Kahan theorem with
$\mathrm{L}^{(p)}=\mathrm{L}^{\infty}+\Delta$ and
$\Delta=\mathrm{L}^{(p)}-\mathrm{L}^{(\infty)}$, for $x_{p}$ and $x_{\infty}$
the eigenvectors of $\mathrm{L}^{(p)}$ and $\mathrm{L}^{(\infty)}$ with
largest eigenvalues respectively,

(29) $\displaystyle|x_{p}-x_{\infty}|\leq\frac{3\hskip
2.0pt}{\delta^{2}}\hskip 2.0pt|\mathrm{L}^{(p)}-\mathrm{L}^{(\infty)}|\hskip
1.0pt.$

Note the condition $\langle x_{p},x_{\infty}\rangle\geq 0$ is without loss of
generality as the orientation of the eigenvectors is always arbitrary. The
following result follows immediately.

######  Proposition 4.3.

Fix $n\geq 2$ and suppose that Assumptions 2.1, 2.2 and 2.3 hold. Then, we
have $|x_{p}-x_{\infty}|\to 0$ as $p\uparrow\infty$ almost surely.

### 4.3. Verifying Assumption 3.5

We make use of Assumptions 2.1–2.2. There are four items to verify for
$\theta=\chi_{n}\beta$ and $w$ in $(\ref{dev})$. Take, $m=\chi_{n}\mu$,
$\xi=\chi_{n}\sigma$ and $\nu=\delta/\sqrt{n}$.

  1. (i)

We have $m_{\infty}(\theta)=\chi_{n}\mu$ and
$\mathrm{s}^{2}_{\infty}(\theta)=\chi_{n}^{2}\hskip 1.0pt\sigma^{2}$ by
Assumption 2.2 part (i).

  2. (ii)

To see that $m_{\infty}(w)=0$, we compute $m(w)$ using $(\ref{dev})$ which
gives

$m(w)=\chi_{n}m(\beta)\langle
x_{p},x_{p}-x_{\infty}\rangle+m(\mathcal{E}x_{p})/\sqrt{n}$

The first term vanishes by Proposition 4.3 since $|\langle
x_{p},x_{p}-x_{\infty}\rangle|\leq|x_{\infty}-x_{p}|$ and $m_{\infty}(\beta)$
and $\chi_{n}$ are finite almost surely by Assumption 2.2 part (ii). The
second term vanishes because $m(\mathcal{E}x_{p})$ may be written as a linear
combination of a fixed $n$ realizations of $m(\epsilon)$ with coefficients
being the entries of $x_{p}$, and $|x_{\infty}|$ is bounded. To this end,
$m(\epsilon)=m_{p}(\epsilon)={\textstyle\sum\nolimits}_{i=1}^{p}\phi_{i}m(\gamma^{i})=m_{p}(\varphi)$
which tends to zero as $p\uparrow\infty$ (almost surely) by Assumption 2.2
part (iii).

Similarly, to verify that $\mathrm{s}^{2}_{\infty}(w)=\nu$, we use
$(\ref{dev})$ again to calculate that

$\mathrm{s}^{2}(w)=\chi_{n}^{2}(\langle\beta,\beta\rangle/p)\langle
x_{p},x_{p}-x_{\infty}\rangle^{2}+\langle\mathcal{E}x_{p},\mathcal{E}x_{p}\rangle/(pn)-m^{2}(w)$

where we have used the fact that $\langle\beta,\mathcal{E}\rangle=0_{n}$ by
Assumption 2.1. Since $\chi_{n}$ and the limit
$\chi_{n}^{2}(\sigma^{2}+\mu^{2})$ of $\langle\beta,\beta\rangle/p$, as above,
by Proposition 4.3 the first term tends to zero as $p\uparrow\infty$. For the
second term, we note that $\langle\mathcal{E}x_{p},\mathcal{E}x_{p}\rangle$
may be written as convex combination of a fixed $n$ realizations of
$\langle\epsilon,\epsilon\rangle$ with coefficients being the entries of
$x_{p}$ squared, and $|x_{p}|^{2}=|x_{\infty}|^{2}=1$. We have
$\langle\epsilon,\epsilon\rangle/p=\sum_{i=1}^{p}\phi_{i}^{2}/p=\mathrm{s}^{2}(\phi)+m^{2}(\phi)$
as the $\\{\gamma^{i}\\}_{i=1}^{p}$ are orthonormal. By Assumption 2.2 parts
(iii) and (iv), we have $\mathrm{s}^{2}_{\infty}(\phi)=\delta^{2}$ and
$m_{\infty}(\phi)=0$. It follows that the second term converges to
$\delta^{2}|x_{\infty}|^{2}/n=\delta^{2}/n=\nu^{2}$. The last term tends to
zero since $m_{\infty}(w)=0$ as above, and the claim now follows.

  3. (iii)

We have $\textup{cov}(\theta,w)=\langle\theta,w\rangle/p-m(\theta)m(w)$ and
since $m_{\infty}(\theta)$ is finite, the second term vanishes as
$m_{\infty}(w)$ as above. Again by $(\ref{dev})$ and since
$\langle\beta,\mathcal{E}\rangle=0_{n}$,

$\displaystyle\langle\theta,w\rangle/p=\chi_{n}^{2}(\langle\beta,\beta\rangle/p)\hskip
1.0pt\langle x_{p},x_{p}-x_{\infty}\rangle$

which vanishes in the limit by the same arguments as in (ii) above.

  4. (iv)

To see that $\hat{\nu}=\hat{\nu}_{p}$ in $(\ref{noisest})$ is an
asymptotically exact estimate of $\nu=\delta/\sqrt{n}$, we use Proposition
4.2. Since
$(\textup{tr}(\mathrm{S})-\mathscr{s}^{2}_{p})/p=\sum_{i=2}^{n}\mathscr{s}^{2}_{p-i+1}/p$,
under the hypotheses of Proposition 4.2 converges almost surely to
$\delta^{2}(n-1)/n$, $\hat{\nu}_{\infty}=\nu$.

##  Appendix A Auxiliary proofs

As shown in Section 4.3, the hypotheses of Proposition 2.5 and Proposition 2.7
(i.e., Assumptions 2.1 –2.4) guarantee that the conditions on
$\\{\theta_{i}\\}_{i\geq 1}$ and $\\{w_{i}\\}_{i\geq 1}$ in Assumption 3.5 are
satisfied. Consequently the proofs of these two results, as well as that of
Proposition 3.7 which requires Assumption 3.5 directly, reduce to the
calculations below. The proof of Lemma 3.6 is omitted as it is elementary and
that of Lemma 3.4 is a direct consequence of some of the expressions below.

For any $\mathrm{c}\in\mathbb{R}$ and
$\eta(\mathrm{c})=m(\eta)+\mathrm{c}\hskip 1.0pt(\eta-m(\eta))$, by direct
calculation

$\displaystyle\textup{MSE}\hskip 1.0pt(\eta(\mathrm{c})\hskip 1.5pt|\hskip
1.5pt\theta)$
$\displaystyle=(m(\eta)-m(\theta))^{2}+\mathrm{s}^{2}(\theta)+\mathrm{c}^{2}\mathrm{s}^{2}(\eta)-2\mathrm{c}\hskip
1.0pt\textup{cov}(\eta,\theta)\hskip 1.0pt.$

When $\mathrm{c}=1$ for which $\eta(1)=\eta$ and applying Lemma 3.6 yields

$\displaystyle\textup{MSE}_{\infty}(\eta\hskip 1.5pt|\hskip 1.5pt\theta)$
$\displaystyle=\xi^{2}+(\xi^{2}+\nu^{2})-2\xi^{2}=\nu^{2}=\delta^{2}/n$

which proves Proposition 2.5. The sine of the angle squared metric is computed
as

$\displaystyle\textup{SPH}\hskip 1.0pt(\eta(\mathrm{c})\hskip 1.5pt|\hskip
1.5pt\theta)=1-\bigg{(}\frac{\langle\eta(\mathrm{c}),\theta\rangle}{|\eta(\mathrm{c})||\theta|}\bigg{)}^{2}=1-\frac{\big{(}m(\eta)m(\theta)+\mathrm{c}\hskip
1.0pt\textup{cov}(\eta,\theta)\big{)}^{2}}{(m^{2}(\eta)+\mathrm{c}^{2}\mathrm{s}^{2}(\eta))(m^{2}(\theta)+\mathrm{s}^{2}(\theta))}\hskip
1.0pt.$

Using the raw estimate $\eta=\eta(1)$ for which $\mathrm{c}=1$ we deduce by
Lemma 3.6 that

$\displaystyle\textup{SPH}_{\infty}(\eta\hskip 1.5pt|\hskip
1.5pt\theta)=1-\frac{m^{2}+\xi^{2}}{m^{2}+\xi^{2}+\nu^{2}}=\frac{\nu^{2}}{\xi^{2}+m^{2}+\nu^{2}}=\frac{\mathrm{r}^{2}_{\infty}}{\textup{SNR}^{2}+\mathrm{r}^{2}_{\infty}}$

which establishes Proposition 2.7 with SNR and $\mathrm{r}_{\infty}$ in
$(\ref{snr})$ (c.f. $(\ref{asnr})$).

Note that minimizing the expressions for $\textup{MSE}\hskip
1.0pt(\eta(\mathrm{c})\hskip 1.5pt|\hskip 1.5pt\theta)$ and
$\textup{SPH}\hskip 1.0pt(\eta(\mathrm{c})\hskip 1.5pt|\hskip 1.5pt\theta)$
above over $\mathrm{c}\in\mathbb{R}$ yields the $\mathrm{c}^{\textup{MSE}}$
and $\mathrm{c}^{\textup{SPH}}$ in $(\ref{cex})$ proving Lemma 3.4. The limits
as $p\uparrow\infty$ of these quantities is easily verified as
$\mathrm{c}_{\infty}=\frac{\textup{SNR}^{2}}{1+\textup{SNR}^{2}}$ in
$(\ref{clim})$ using Lemma 3.6. This establishes the last part of Proposition
3.7. To prove the first part, we again apply Lemma 3.6 to deduce that
$\mathrm{c}=\mathrm{c}_{p}=1-\frac{\hat{\nu}_{p}}{\mathrm{s}^{2}_{p}(\eta)}\sim
1-\frac{\nu^{2}}{\xi^{2}+\nu^{2}}=\mathrm{c}_{\infty}$ as above. Also,

$\displaystyle\textup{MSE}_{\infty}(\eta(\mathrm{c})\hskip 1.5pt|\hskip
1.5pt\theta)$
$\displaystyle=\xi^{2}+\mathrm{c}^{2}_{\infty}(\xi^{2}+\nu^{2})-2\mathrm{c}_{\infty}\xi^{2}$
$\displaystyle=\xi^{2}+\mathrm{c}_{\infty}\xi^{2}-2\mathrm{c}_{\infty}\xi^{2}$
$\displaystyle=\xi^{2}(1-\mathrm{c}_{\infty})$
$\displaystyle=\frac{\xi^{2}}{1+\textup{SNR}^{2}}=\mathrm{c}_{\infty}\nu^{2}=\mathrm{c}_{\infty}\textup{MSE}_{\infty}(\eta\hskip
1.5pt|\hskip 1.5pt\theta)$

as required. Similarly, by Lemma 3.6, we derive that

$\displaystyle\textup{SPH}_{\infty}(\eta(\mathrm{c})\hskip 1.5pt|\hskip
1.5pt\theta)$
$\displaystyle=1-\frac{(m^{2}+\xi^{2}\mathrm{c}_{\infty})^{2}}{(m^{2}+\mathrm{c}_{\infty}^{2}(\xi^{2}+\nu^{2}))(\mu^{2}+\sigma^{2})}=(1-\mathrm{c}_{\infty})\hskip
1.0pt\mathrm{r}^{2}_{\infty}$
$\displaystyle=(1-\mathrm{c}_{\infty})(\textup{SNR}^{2}+\mathrm{r}^{2}_{\infty})\textup{SPH}_{\infty}(\eta\hskip
1.5pt|\hskip 1.5pt\theta)$
$\displaystyle=\bigg{(}\frac{\textup{SNR}^{2}+\mathrm{r}^{2}_{\infty}}{1+\textup{SNR}^{2}}\bigg{)}\textup{SPH}_{\infty}(\eta\hskip
1.5pt|\hskip 1.5pt\theta)=\mathrm{d}_{\infty}\textup{SPH}_{\infty}(\eta\hskip
1.5pt|\hskip 1.5pt\theta)$

for $\mathrm{d}_{\infty}$ as in $(\ref{slim})$. This concludes the proof of
Proposition 3.7.

## References

  * [1] [Ahn et al.]Ahn, Marron, Muller Chi2007ahn2007 Ahn, J., Marron, J., Muller, K. M. Chi, Y.-Y. 2007, ‘The high-dimension, low-sample-size geometric representation holds under mild conditions’, Biometrika 94(3), 760–766.
  * [2] Anderson2003anderson2003 Anderson, T. W. 2003, An Introduction to Multivariate Statistical Analysis, John Wiley & Sons, Inc.
  * [3] [Aoshima et al.]Aoshima, Shen, Shen, Yata, Zhou Marron2018aoshima2018 Aoshima, M., Shen, D., Shen, H., Yata, K., Zhou, Y.-H. Marron, J. 2018, ‘A survey of high dimension low sample size asymptotics’, Australian & New Zealand journal of statistics 60(1), 4–19.
  * [4] Bock1975bock1975 Bock, M. E. 1975, ‘Minimax estimators of the mean of a multivariate normal distribution’, The Annals of Statistics pp. 209–218.
  * [5] Brown1971brown1971 Brown, L. D. 1971, ‘Admissible estimators, recurrent diffusions, and insoluble boundary value problems’, The Annals of Mathematical Statistics 42(3), 855–903.
  * [6] Brown Zhao2012brown2012 Brown, L. D. Zhao, L. H. 2012, ‘A geometrical explanation of Stein shrinkage’, Statistical Science pp. 24–30.
  * [7] [Donoho et al.]Donoho, Gavish Johnstone2018donoho2018 Donoho, D. L., Gavish, M. Johnstone, I. M. 2018, ‘Optimal shrinkage of eigenvalues in the spiked covariance model’, Annals of statistics 46(4), 1742.
  * [8] Efron Morris1975efron1975 Efron, B. Morris, C. 1975, ‘Data analysis using stein’s estimator and its generalizations’, Journal of the American Statistical Association 70(350), 311–319.
  * [9] [Fan et al.]Fan, Guo Zheng2020fan2020 Fan, J., Guo, J. Zheng, S. 2020, ‘Estimating number of factors by adjusted eigenvalues thresholding’, Journal of the American Statistical Association pp. 1–10.
  * [10] [Fourdrinier et al.]Fourdrinier, Strawderman Wells2018fourdrinier2018 Fourdrinier, D., Strawderman, W. E. Wells, M. T. 2018, Shrinkage estimation, Springer.
  * [11] [Goldberg et al.]Goldberg, Papanicalaou Shkolnik2021goldberg2021 Goldberg, L., Papanicalaou, A. Shkolnik, A. 2021, Dispersion bias.  Working paper.
  * [12] [Goldberg et al.]Goldberg, Papanicolaou, Shkolnik Ulucam2020goldberg2020 Goldberg, L. R., Papanicolaou, A., Shkolnik, A. Ulucam, S. 2020, ‘Better betas’, The Journal of Portfolio Management 47(1), 119–136.
  * [13] Gruber2017gruber2017 Gruber, M. H. 2017, Improving efficiency by shrinkage: the James–Stein and ridge regression estimators, Routledge.
  * [14] Gupta Peña1991gupta1991 Gupta, A. K. Peña, E. A. 1991, ‘A simple motivation for James-Stein estimators’, Statistics & probability letters 12(4), 337–340.
  * [15] [Hall et al.]Hall, Marron Neeman2005hall2005 Hall, P., Marron, J. S. Neeman, A. 2005, ‘Geometric representation of high dimension, low sample size data’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67(3), 427–444.
  * [16] Hoerl Kennard1970hoerl1970 Hoerl, A. E. Kennard, R. W. 1970, ‘Ridge regression: Biased estimation for nonorthogonal problems’, Technometrics 12(1), 55–67.
  * [17] Horn Johnson2013horn2012 Horn, R. A. Johnson, C. R. 2013, Matrix analysis, Cambridge university press.
  * [18] James Stein1961james1961 James, W. Stein, C. 1961, Estimation with quadratic loss, in ‘Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics’, The Regents of the University of California.
  * [19] Jolliffe Cadima2016jolliffe2016 Jolliffe, I. T. Cadima, J. 2016, ‘Principal Component Analysis: a review and recent developments’, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374(2065), 20150202.
  * [20] Jung Marron2009jung2009 Jung, S. Marron, J. S. 2009, ‘PCA consistency in high dimension, low sample size context’, The Annals of Statistics 37(6B), 4104–4130.
  * [21] [Jung et al.]Jung, Sen Marron2012jung2012 Jung, S., Sen, A. Marron, J. 2012, ‘Boundary behavior in high dimension, low sample size asymptotics of PCA’, Journal of Multivariate Analysis 109, 190–203.
  * [22] Ledoit Wolf2004ledoit2004 Ledoit, O. Wolf, M. 2004, ‘A well-conditioned estimator for large-dimensional covariance matrices’, Journal of multivariate analysis 88(2), 365–411.
  * [23] Paul2007paul2007 Paul, D. 2007, ‘Asymptotics of sample eigenstructure for a large dimensional spiked covariance model’, Statistica Sinica pp. 1617–1642.
  * [24] [Shen et al.]Shen, Shen, Zhu Marron2016shen2016 Shen, D., Shen, H., Zhu, H. Marron, J. 2016, ‘The statistics and mathematics of high dimension low sample size asymptotics’, Statistica Sinica 26(4), 1747.
  * [25] Stein1956stein1956 Stein, C. 1956, Inadmissibility of the usual estimator for the mean of a multivariate normal distribution, in ‘Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics’, The Regents of the University of California.
  * [26] Stein1986stein1986 Stein, C. 1986, ‘Lectures on the theory of estimation of many parameters’, Journal of Soviet Mathematics 34(1), 1373–1403.
  * [27] Stigler1990stigler1990 Stigler, S. M. 1990, ‘The 1988 Neyman memorial lecture: a Galtonian perspective on shrinkage estimators’, Statistical Science pp. 147–155.
  * [28] Tibshirani1996tibshirani1996 Tibshirani, R. 1996, ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society: Series B (Methodological) 58(1), 267–288.
  * [29] Wang Fan2017wang2017 Wang, W. Fan, J. 2017, ‘Asymptotics of empirical eigenstructure for high dimensional spiked covariance’, Annals of statistics 45(3), 1342.
  * [30] Weyl1912weyl1912 Weyl, H. 1912, ‘Das asymptotische verteilungsgesetz der eigenwerte linearer partieller differentialgleichungen (mit einer anwendung auf die theorie der hohlraumstrahlung)’, Mathematische Annalen 71(4), 441–479.
  * [31] [Yu et al.]Yu, Wang Samworth2015yu2015 Yu, Y., Wang, T. Samworth, R. J. 2015, ‘A useful variant of the Davis–Kahan theorem for statisticians’, Biometrika 102(2), 315–323.
  * [32] Zou Hastie2005zou2005 Zou, H. Hastie, T. 2005, ‘Regularization and variable selection via the elastic net’, Journal of the royal statistical society: series B (statistical methodology) 67(2), 301–320.
  * [33]