# CHARGE-INJECTION DEVICE IMAGING OF SIRIUS WITH CONTRAST RATIOS
GREATER THAN 1:26 MILLION

Sailee M. Sawant Department of Aerospace, Physics and Space Sciences, 150 W.
University Blvd., Florida Institute of Technology,
Melbourne, FL 32901, USA  Daniel Batcheldor Department of Aerospace, Physics
and Space Sciences, 150 W. University Blvd., Florida Institute of Technology,
Melbourne, FL 32901, USA  Southeastern Universities Research Association,
Laboratory Support Services and Operations,
Kennedy Space Center, FL 32899, USA

###### Abstract

The intrinsic nature of many astronomical objects, such as binary systems,
exoplanets, circumstellar and debris disks, and quasar host galaxies,
introduces challenging requirements for observational instrumentation and
techniques. In each case, we encounter situations where the light from bright
sources hampers our ability to detect surrounding fainter targets. To explore
all features of such astronomical scenes, we must perform observations at the
maximum possible contrast ratios. Charge-injection devices (CIDs) are capable
of potentially exceeding contrast ratios of $\log_{10}{(CR)}>9$ (i.e., 1 part
in 1 billion) due to their unique readout architectures and inherent anti-
blooming abilities. The on-sky testing of the latest generation of CIDs, the
SpectraCAM XDR, has previously demonstrated direct contrast ratios in excess
of 1 part in 20 million from sub-optimal ground-based astronomical
observations that imposed practical limits on the maximum achievable contrast
ratios. Here, we demonstrate the extreme contrast ratio imaging capabilities
of the SXDR using observations of Sirius with the 1.0-m Jacobus Kapteyn
Telescope, La Palma, Spain. Based on wavelet-based analysis and precise
photometric and astrometric calibrations, we report a direct contrast ratio of
$\Delta m_{r}=18.54$, $\log_{10}{(CR)}=7.41\pm 0.08$, or $1$ part in $26$
million. This shows a $29\%$ increase in the achievable contrast ratios
compared to the previous results.

Astronomy Data Analysis, Astronomical Instrumentation, Astronomical
Techniques, Direct Imaging

##  1 Introduction

Direct imaging of fainter targets in the vicinity of bright sources imposes
limitations on the type of contrast observations achievable with ground- and
space-based telescopes. The combined effect of small angular separations and
large magnitude differences severely impacts the possibility of direct
detection of fainter targets (Marois et al., 2008; Oppenheimer & Hinkley,
2009). Additionally, the light from bright sources saturates conventional
imaging instrumentations, e.g., charge-coupled devices (CCDs) and
complementary metal-oxide-semiconductor (CMOS) detectors. Due to their limited
full well depth and dynamic range, the directly achievable contrast ratios of
these 16-bit detectors are restricted to $\log_{10}{(CR)}<5$ (i.e., 1 part in
100 thousand).

Several point-spread function (PSF) suppression techniques are implemented to
mitigate bright source signals and subsequently achieve higher contrast ratios
(i.e., $5<\log_{10}{(CR)}<7$). Currently employed techniques, such as
coronagraphy (e.g., Schneider et al., 2001), nulling interferometry (e.g.,
Bracewell & MacPhie, 1979; Linfield, 2003), integral field spectral
deconvolution (e.g., Sparks & Ford, 2002), ground-based angular differential
imaging (e.g., Marois et al., 2006), and space-based roll subtraction (e.g.,
Lowrance et al., 2005; Schneider et al., 2010), have complex operational
requirements that make faint signal detection a challenge. Additionally, these
techniques are optimized when stable PSFs are achieved.

An alternative is to use appropriate instrumentation that can carry out direct
extreme contrast ratio (ECR) imaging free from the limitations of currently
employed techniques. Charge-injection devices (CIDs) have the intrinsic
ability to achieve extreme contrast ratios (ECRs) of $\log_{10}{(CR)}>9$
(i.e., 1 part in 1 billion) owing to their unique readout architectures and
inherent anti-blooming abilities (Bhaskaran et al., 2008). A preliminary study
of the Sirius field using the latest generation of commercially available
CIDs, the SpectraCAM XDR (SXDR), at the Florida Tech 0.8-m Ortega telescope
demonstrated a direct raw contrast ratio of $\Delta m_{v}=18.3$,
$\log{(CR)}=7.3$, or 1 part in 20 million at an angular separation of two
arcminutes (Batcheldor et al., 2016). However, the atmospheric conditions in
Florida imposed practical limitations on the type of contrast ratios that are
otherwise achievable using CIDs. Nonetheless, this result provided a
motivation to carry out a direct CID-based ECR imaging of Sirius using the
1.0-m Jacobus Kapteyn Telescope (JKT) at the Roque de Los Muchachos
Observatory, La Palma, Spain.

In this paper, we report the results of the direct ECR imaging of Sirius using
the SXDR provided by ThermoFisher Scientific Inc. In $\S$2, we provide a brief
introduction to CIDs and discuss the details of our observations. $\S$3
describes the implemented data reduction and analysis. In $\S$4, we present
the results of the direct ECR imaging of Sirius. In $\S$5, we discuss these
results and outline the future prospects for implementing CIDs as astronomical
imaging detectors. $\S$6 concludes.

##  2 Instrumentations and Observations

CIDs are an array of X-Y addressable photosensitive metal-oxide-silicon (MOS)
capacitors (Michon, G., 1974). First developed in 1970, the original CID
imagers incorporated complex on-chip readout electronics and external (i.e.,
off-chip) amplifiers that delivered a high read noise of $\sim 600\,e^{-}$
(McCreight & Goebel, 1981; McCreight et al., 1986). As a result, CCDs (read
noise of $<10\,e^{-}$) have been preferred as astronomical detectors. However,
modern CIDs are fabricated with pre-amplifier per pixel (PPP) architectures
(active pixel technology), random access decoders, and non-destructive readout
(NDRO) mechanisms to deliver read noise levels comparable to CCDs (Ninkov et
al., 1994; Eid, 1995; Kimble et al., 1995; Bhaskaran et al., 2008).

Modern CIDs are also designed to support random access integration (RAI) with
extreme dynamic range (XDR). The RAI mode utilizes a random access integrator
to identify, monitor, and control the data acquisition on intensely
illuminated pixels. A region of interest (ROI) is identified in a short pre-
exposure (e.g., 0.1 second). The pixels within the ROI are integrated until an
NDRO signal threshold (i.e., $75\%$ of the full well capacity) is reached. The
detected signal is readout and stored with a timestamp. The accumulated charge
is injected into the underlying substrate, and the pixel operation is
repeated. The data is then combined to create a final image with a wider
dynamic range. Since the values of intensely illuminated pixels do not
restrict the full well depths of CIDs, the pixels may be integrated, injected,
and re-integrated multiple times during a pre-determined exposure. This allows
CIDs to extend the upper limit of the dynamic range, generate 32-bit images,
and deliver a potential contrast ratio of $\Delta m_{v}\sim 24$,
$\log_{10}{(CR)}\sim 9.6$, or 1 part in 4 billion (Bhaskaran et al., 2008).

Due to their unique pixel architecture, true random accessibility, and NDRO
capabilities, CIDs offer several advantages over CCDs: negligible charge-
transfer losses, extreme radiation tolerance ($>5\,$Mrad), large full well
capacity, high fill-factor, broad spectral response (156 nm - 1100 nm), and
intrinsic anti-blooming abilities (Bhaskaran et al., 2008). Therefore, CIDs
are excellent prospects to perform direct ECR imaging with ground- and space-
based telescopes.

The data presented here was acquired using the SXDR installed on the 1.0-m
JKT, a facility operated by the Southeastern Association for Research in
Astronomy (SARA) at the Roque de Los Muchachos Observatory, La Palma, Spain.
The SXDR utilizes a $2048\times 2048$ $12\,\mu$m pixel Random Access CID820
(RACID) sensor with PPP architecture. The field of view is $10^{\prime}\times
10^{\prime}$ with an image scale of $0\farcs 30$/pixel. The pixels have a
linear (within $2\%$) full well capacity of $268,000\,e^{-}$. They saturate at
$305,000\,e^{-}$. The peak quantum efficiency at $525\,$nm ($\sim$ V-band) is
$48\%$ (Bhaskaran et al., 2008). The detector is hermetically sealed and
cooled to $-45.6^{\circ}$C using an ethylene glycol re-circulation system. At
this temperature, the dark current is $5\,e^{-}/\text{s}$ with a conversion
gain of $6.2\,e^{-}/$ADU. The noise associated with a single read is
$44\,e^{-}$ RMS. This reduces to $5.8\,e^{-}$ RMS with 128 NDROs for the
sampling frequency of 2.1 MHz (Batcheldor et al., 2016).

We conducted observations of Sirius on March 8th, 2017, at an airmass of 1.4.
We also observed the Landolt standard star field of SA98 670. All images are
acquired with the Johnson-Cousins B, V, R, and I filters. The average DIMM
(Differential Image Motion Monitor) seeing was $1\farcs 7$ FWHM. The moon
illumination was $86\%$.

##  3 Data Reduction and Analysis

Figure 1: Pre-processed I-band SXDR image of Sirius with an exposure time of
180 seconds. The field of view is $10^{\prime}\times 10^{\prime}$. The signal
from Sirius (RA = $06^{\text{h}}45^{\text{m}}08.26^{\text{s}}$, DEC =
$-16^{\circ}43^{\prime}19\farcs 03$) is not saturated. The standard pre-
processing routines removed instrumental imperfections; however, it did not
suppress the visible artifacts, such as the multiple annular patterns,
atmospheric halo, and diffraction spikes.

Figure 1 shows a pre-processed I-band SXDR image of the $10^{\prime}\times
10^{\prime}$ field around Sirius at an exposure time of 180 seconds. The
signal from Sirius (RA = $06^{\text{h}}45^{\text{m}}08.26^{\text{s}}$, DEC =
$-16^{\circ}43^{\prime}19\farcs 03$) is not saturated. The pre-processing
routines successfully removed instrumental imperfections due to dark and flat-
field responses; however, an additional post-processing routine was required
to suppress the visible artifacts. The internal reflections from the
apochromatic field corrector lenses in the JKT produced multiple ghosts (i.e.,
annular patterns). It is worth noting that such ghosts are absent in the
preliminary observations of Sirius (see Figure 4 in Batcheldor et al., 2016).
These artifacts are purely intrinsic to the JKT due to its optical design to
acquire a flatter focal plane and a wider field of view (Keel et al., 2016).
In addition to this, the atmospheric halo and diffraction spikes dominated the
sky background. These effects imposed limitations on the full dynamic range of
the image. Therefore, we implemented an additional post-processing routine to
mitigate the effects of ghosts, atmospheric halo, and diffraction spikes.

###  3.1 Wavelet Transformation

For the cases where the background level varies uniformly, a simple median
filter is applied to generate a background model, which is then subtracted
from pre-processed images (e.g., Bertin & Arnouts, 1996; Bradley et al.,
2020). Additionally, appropriate sigma-clipping techniques are used to achieve
better accuracy (Starck & Murtagh, 2002a). Diffused halos from bright sources
are suppressed using several methods, such as radial-gradient methods (e.g.,
Bonnet-Bidaud et al., 2000) and isophotal analysis (e.g., Batcheldor et al.,
2016).

However, the presence of multiple ghosts in our data (as seen in Figure 1)
resulted in a non-uniform background with abrupt spatial variations. The
amalgamation of these patterns with faint source signals highlighted an issue
of complex hierarchical structures. Therefore, we applied the àtrous wavelet
transform to decompose and analyze the SXDR images at different resolution-
related scales (Starck & Murtagh, 2002b).

The àtrous wavelet transform yields a discrete, stationary, isotropic, and
shift-invariant transformation, which makes it well-suited for astronomical
image and data analysis (Starck & Murtagh, 2002b; Mertens & Lobanov, 2015). We
followed the algorithm discussed in Starck & Murtagh (2002a) and applied the
àtrous wavelet transform to the SXDR images. This yielded a set of resolution-
related views for a given input SXDR image, which is described as follows:

$f(x,y)=\sum_{j=1}^{J}w_{j}(x,y)+c_{J}(x,y),$ (1)

where $f(x,y)$ is a mathematical representation of the input SXDR image,
$w_{j}(x,y)$ is a set of resolution-related views of the image, called wavelet
scales, and $c_{J}$ is the last smoothed array. The set
$w=\mathcal{W}f(x,y)=\\{w_{1}(x,y),...,w_{J}(x,y),c_{J}(x,y)\\}$ is known as
the wavelet transform of the input SXDR image (Mertens & Lobanov, 2015).

We modeled the noise in the wavelet space to accurately assess the background
variation across the input SXDR image. At each resolution scale $j$, we
obtained a wavelet coefficient distribution and carried out a statistical
significance test to determine whether a given wavelet coefficient
$w_{j}(x,y)$ is due to signal (i.e., significant) or noise (i.e., not
significant). We defined our null hypothesis, $\mathcal{H}_{0}$, as a value
that is locally constant across the input SXDR image. We rejected the
hypothesis $\mathcal{H}_{0}$ when the corresponding $p$-value is less than a
pre-specified significance level, $\alpha$ (Miller et al., 2001). Several
thresholding techniques exist for multiple hypothesis testing (e.g.,
$2\sigma$, $3\sigma$, and the Bonferoni methods). However, these techniques
fail to control a critical quantity, the false discovery rate (FDR), the
fraction of false rejections based on the total number of rejections (Miller
et al., 2001). The FDR method outmatches the standard techniques by
successfully controlling this quantity (Hopkins et al., 2002). The
corresponding FDR ratio is defined as follows:

$\displaystyle\textrm{FDR}=\frac{N^{\mathrm{reject}}_{\mathrm{null\;true}}}{N^{\mathrm{reject}}}$
(2)

where $N^{\mathrm{reject}}_{\mathrm{null\;true}}$ is the number of pixels that
are truly not significant but declared significant (i.e., false positives) and
$N^{\mathrm{reject}}$ is the total number of pixels declared significant. This
method ensures that, on average, the FDR is always less than or equal to
$\alpha$ (Miller et al., 2001).

We selected an $\alpha$ value such that $0\leq\alpha\leq 1$. At each scale, we
calculated the $p$-values of the $N$ coefficients, sorted from smallest to
largest. We defined and evaluated a quantity, $d$, such that:

$\displaystyle
d=max\left\\{i:P_{i}<P_{D}\right\\},\mathrm{\;where\;}P_{D}=\frac{i\alpha}{c_{N}N}$
(3)

where $P_{i}$ is the $p$-value of the ith ordered coefficients and $c_{N}$ is
a constant, which is equal to 1 when the $p$-values are statistically
independent. We rejected all hypotheses whose $p$-values are less or equal to
$P_{D}$. We then selected the FDR threshold as the wavelet coefficient
$w(x,y)$ with the largest index $i$ whose $p$-value satisfied the above
inequality. In this way, we calculated an adaptive significant threshold at
each scale. A detailed description of this procedure is discussed in Miller et
al. (2001).

We followed the reconstruction algorithm presented in Starck & Murtagh (2002a)
and produced a filtered version of the input SXDR image. To avoid the loss of
flux due to thresholding, we improved our filtering algorithm by using the Van
Cittert iteration algorithm (Starck et al., 1998). We determined the wavelet
transform of the filtered SXDR image and ensured that it produced the same
significant wavelet coefficients. For each iteration, we calculated the
residual signal. We determined the wavelet transform of the residual signal,
applied thresholding, and reconstructed the threshold error signal. The
residual signal was added back to the filtered SXDR image. We defined our
convergence criteria by calculating the estimated relative error and comparing
it to a pre-specified tolerance value, $\epsilon=0.002$.

###  3.2 Faint Source Detection

We used the detection program photutils.detection, an affiliated package of
AstroPy (Astropy Collaboration et al., 2018; Bradley et al., 2020), to
identify sources in the reconstructed noise-reduced SXDR images. This program
uses the DAOFIND algorithm (Stetson, 1987) to identify local maxima against a
user-specified threshold value. To detect the faintest sources, we applied the
same detection algorithm to the first filtered resolution-view of the SXDR
images. This is only possible in wavelet space because the large-scale
features, such as the visible artifacts seen in Figure 1, are suppressed at
the smallest scale revealing the faint small-scale sources (Starck & Murtagh,
2002c). We extracted the centroid coordinates (x, y) of the detected sources
with sub-pixel accuracy. Furthermore, we constrained our detection criteria to
extract sources that have an ellipticity of 1.0 (i.e., point sources).

###  3.3 Photometric and Astrometric Calibrations

The process of transforming our instrumental magnitudes to the standard UBVRI
magnitudes was carried out using the Gaia Early Data Release 3 (Gaia EDR3)
catalog data (Gaia Collaboration et al., 2016; Riello et al., 2021). We found
our photometric accuracy using the B-, V-, R-, and I-band images of the
Landolt standard star field of SA98 670. The resulting mean absolute
percentage errors (MAPEs) in the transformed $B,V,R,$ and $I$ magnitudes of
the standard stars of SA98 670 are $0.167\%,0.048\%,0.072\%,\textrm{ and
}0.068\%$, respectively.

To match the positions of our detected sources with previously known sources,
we acquired the astrometric coordinates of the Gaia EDR3 cataloged sources
within the same field of view. We cross-correlated pixel coordinates of the
Gaia EDR3 sources with centroids of the detected sources. We compared the
cataloged coordinates of the cross-correlated sources with the corresponding
calculated coordinates. The resulting mean absolute accuracies in the right
ascension (RA) and declination (DEC) are $0\farcs 321$ and $0\farcs 258$,
respectively.

##  4 Results

Figure 2: Post-processed I-Band SXDR image of the Sirius field with an
exposure time of 180 seconds. The SIMBAD and the detected Gaia EDR3 cataloged
sources are marked in red and blue, respectively. Figure 3: Contrast ratios of
the detected sources in the B-, V-, R-, and I-band as a function of angular
separation (in arcseconds) from Sirius. Detected sources with SNRs $<10$,
$10\leq\textrm{SNRs}<20$, and SNRs $\geq 20$ are shown in red, blue, and
green, respectively. For sources with SNRs $<10$, we report their raw contrast
ratios. Smaller errors are contained within individual markers.

Filter |  Source ID |  Gaia |  |  |
---|---|---|---|---|---
|  Calibrated Magnitude (mag) |  $\log_{10}(CR)$ |  SNR |  |
B | 1 | 2947057338479481600 | 15.216 $\pm$ 0.079 | 6.67 $\pm$ 0.03 | 48.66
V | 2 | 2947046961838568064 | 16.307 $\pm$ 0.189 | 7.11 $\pm$ 0.08 | 15.53
| 3 | 2947046686960684288 | 15.587 $\pm$ 0.142 | 6.82 $\pm$ 0.06 | 32.46
R | 4 | 2947063042187126400 | 17.077 $\pm$ 0.203 | 7.41 $\pm$ 0.08 | 7.39
| 5 | 2947062939108083456 | 16.753 $\pm$ 0.166 | 7.29 $\pm$ 0.07 | 11.54
| 6 | 2947051188086163584 | 16.115 $\pm$ 0.113 | 7.03 $\pm$ 0.05 | 24.98
I | 7 | 2947063759445801472 | 16.210 $\pm$ 0.170 | 7.06 $\pm$ 0.07 | 7.73
| 8 | 2947062973476413824 | 15.560 $\pm$ 0.102 | 6.80 $\pm$ 0.04 | 12.36
| 9 | 2947049882416355712 | 14.837 $\pm$ 0.063 | 6.51 $\pm$ 0.03 | 36.66
Table 1: Some of the detected faint Gaia EDR3 cataloged sources in the Sirius
field. These sources demonstrate contrast ratios at different SNR levels in
each band. No sources with SNRs $<10\textrm{ and}<20$ are detected in the V-
and B-band, respectively.

Figure 2 shows the post-processed I-band SXDR image of Sirius. The SIMBAD and
detected Gaia EDR3 cataloged sources are marked in red and blue, respectively.

We calculated direct contrast ratios based on the differences between the
calibrated apparent magnitudes of the detected sources and the cataloged
magnitudes of Sirius (i.e., $\Delta m$). These contrast ratios are usually
reported as follows:

$\log_{10}(CR)=-\frac{\Delta m}{2.5}$ (4)

Figure 3 shows direct contrast ratios of the detected Gaia EDR3 cataloged
sources as a function of angular separation (in arcseconds) from Sirius. Since
the Gaia EDR3 catalog does not include data on Sirius, all contrast ratios and
angular separations are calculated with respect to the SIMBAD cataloged data
of Sirius (Wenger et al., 2000). We acquired the maximum direct contrast ratio
of $\Delta m_{r}=18.54$, $\log_{10}{(CR)}=7.41\pm 0.08$, or $1$ part in $26$
million, a $29\%$ increase from the raw contrast ratio of
$\log_{10}{(CR)}=7.3$ demonstrated from Florida (Batcheldor et al., 2016). For
the cases where the Gaia EDR3 cataloged sources are detected with SNRs $<10$,
we calculated their raw contrast ratios (i.e., contrast ratios based on their
cataloged magnitudes). These contrast ratios are marked in red.

Table 4 lists some of the detected faint Gaia EDR3 cataloged sources. These
sources are detected at different confidence levels based on their SNRs. The
corresponding $8^{\prime\prime}\times 8^{\prime\prime}$ thumbnails are shown
in Figure 4. No sources with SNRs $<10$ and $<20$ are detected in the V- and
B-band, respectively.

Figure 4: $8^{\prime\prime}\times 8^{\prime\prime}$ thumbnails for the
detected faint Gaia EDR3 cataloged sources listed in Table 4.

##  5 Discussion

Due to the extreme brightness of Sirius, previous ground- and space-based
observations have prompted the need for appropriate imaging instrumentation
and techniques. For instance, Bonnet-Bidaud & Gry (1991) and Bonnet-Bidaud et
al. (2000) used specially designed coronagraphic devices on ground-based
telescopes to look for companions around Sirius. Bond et al. (2017) determined
several orbital parameters of the Sirius system using an extensive volume of
Hubble Space Telescope (HST) observations and ground-based photographic
images. Even with HST, significant observing challenges were encountered; no
combination of a filter and short exposure time resulted in an unsaturated
image of Sirius (Bond et al., 2017). However, as we have demonstrated here, it
is now possible to acquire an unsaturated signal from Sirius and detect
sources at comparable inner working angles (IWAs) using a 1.0-m telescope at
an exposure time of 180 seconds without imposing complex operational
requirements.

We also detected Sirius B at an angular separation of $8\farcs 401$ from
Sirius. Figure 5 shows a zoomed-in view of the post-processed I-Band SXDR
image of the field around Sirius. The annotated source is Sirius B. The
separation from Sirius is consistent with the proper motion data (e.g., Wenger
et al., 2000; Gaia Collaboration et al., 2016). This is an improvement in the
IWA of the SXDR because Sirius B, despite being $8^{\textrm{th}}$ magnitude,
was not detected from Florida (Batcheldor et al., 2016).

Figure 5: Central region of the post-processed I-band SXDR image of the Sirius
field. Sirius B, annotated in red, is at an angular separation of $8\farcs
401$ from Sirius. This is consistent with the proper motion data (e.g., Wenger
et al., 2000; Gaia Collaboration et al., 2016).

Within the $5^{\prime}$ radial field around Sirius, the SIMBAD catalog reports
a total of 47 sources out of which [BG91] 1-9 (Bonnet-Bidaud & Gry, 1991) are
duplicated with [BCL2000] 1-9 (Bonnet-Bidaud et al., 2000). This is because
Bonnet-Bidaud & Gry (1991) defined spatial positions as offsets from Sirius,
whereas Bonnet-Bidaud et al. (2000) reported the absolute positions.
Furthermore, there is a 15-year gap between these observations. We applied
corrections to address these source duplicities. Upon analyzing the positions
of these sources over the given baseline, we found out that these sources
exhibit negligible proper motion. However, the proper motion of Sirius
introduced additional systematic offsets to the spatial positions of [BG91].
In addition to this, the photometric and astrometric uncertainties in the
[BG91] data are $\pm 0.3$ mag and $\pm 0\farcs 7$, respectively (Bonnet-Bidaud
& Gry, 1991).

Recently, the Gaia EDR3 catalog reports a total of $1181$ sources (without
duplicities) within the central $10^{\prime}$ square field around Sirius.
These sources have all-inclusive Gaia photometric and astrometric data with
higher precision (see Table 3 in Gaia Collaboration et al. (2021) for
statistical details on uncertainties). Therefore, we used the Gaia EDR3
cataloged sources for our photometric and astrometric calibrations.

To test the accuracy of our analysis, we applied our reduction methods on a
simulated bright star field image. This simulated field consisted of a bright
central source along with 99 point sources. We calculated the contrast ratios
of these simulated sources with respect to the magnitude of the simulated
bright source in a noise-free image. We then added multiple ghost patterns
across the image to replicate the non-uniform background. Lastly, we applied
our reduction methods and recovered the contrast ratios of the simulated point
sources.

Figure 6: Plot of $\log_{10}(CR)$ for a set of simulated sources. The contrast
ratios and angular separations (in arcseconds) are calculated with respect to
a simulated bright source. The unfilled square markers are contrast ratios of
pre-defined simulated sources in a noise-free image. The filled circular
markers are recovered contrast ratios of simulated sources from a noisy image.
The mean error in the recovered contrast ratios is $\pm 0.004$. The horizontal
dashed colored lines show the sky-limited contrast ratios of
$\log_{10}{(CR)}=8.50\pm 0.02,8.63\pm 0.01,8.62\pm 0.01,$ and $8.01\pm 0.01$
in the B-band (marked in blue), V-band (marked in green), R-band (marked in
red) and I-band (marked in purple), respectively.

Figure 6 shows a plot of contrast ratios of simulated sources as a function of
angular distance (in arcseconds) from the simulated bright source. The
unfilled square markers are contrast ratios of pre-defined simulated sources
in a noise-free image. The filled circular markers are recovered contrast
ratios of these simulated sources from a noisy image. The mean error in the
recovered contrast ratios is $\pm 0.004$. The acquired MAPE between the
simulated and recovered contrast ratios is $0.016\%$.

Although we successfully recovered simulated contrast ratios greater than
$\log_{10}{(CR)}=8.0$ (see Figure 6), we did not acquire comparable contrast
ratios from our Sirius observations. Several instrumental and observational
factors imposed practical limitations on the maximum achievable contrast
ratios. The optical design of the JKT produced multiple ghost patterns. The
moon illumination at the time of data collection increased the overall sky
brightness. The resulting non-uniform background and bright sky introduced the
limiting magnitudes of $19.793\pm 0.040,20.120\pm 0.028,20.082\pm
0.018,\textrm{ and }18.818\pm 0.017$ in the B-, V-, R-, and I-band,
respectively. The horizontal dashed lines in Figure 6 represent the
corresponding sky-limited contrast ratios of $\log_{10}{(CR)}=8.50\pm
0.02,8.63\pm 0.01,8.62\pm 0.01,$ and $8.01\pm 0.01$ in the B-band (marked in
blue), V-band (marked in green), R-band (marked in red), and I-band (marked in
purple), respectively. As clearly seen in Figure 3, we observed fewer sources
at shorter wavelengths. This is consistent with the stellar mass function,
which is dominated by the latter spectral types. Furthermore, no sources with
$\textrm{SNRs}<10$ and $<20$ are detected in the V- and B-band, respectively.
Since Sirius is an A1V type star, it is brighter through shorter wavelength
filters. Consequently, Rayleigh scattering increased the standard deviation in
the background noise level, making it challenging to detect fainter sources.

Nevertheless, we identified several sources with SNRs lower than the minimum
threshold. As a result, although visible, these sources did not satisfy our
detection requirements. For instance, in the R-band, we identified Gaia EDR3
2947063454512704768, 2947049985495334144, and 2947056616924962688 that have
direct raw contrast ratios of $\log_{10}{(CR)}=7.68\pm 0.01,7.72\pm 0.01,$ and
$7.73\pm 0.01$, respectively. Under ideal seeing conditions, it would be
possible to detect such sources at $2\sigma$ or $3\sigma$ above the background
level.

CIDs are capable of exceeding contrast ratios of 1 part in 1 billion
(Bhaskaran et al., 2008). Even in the presence of practical limitations, we
achieved a photometric precision of $\pm 0.019,\pm 0.049,\pm 0.041,$ and $\pm
0.044$ mag for $B,V,R,$ and $I$ magnitudes, respectively. As discussed
earlier, the contrast ratios acquired from the simulated bright star field
data show that the implemented reduction and analysis methods are accurate
within a MAPE of $0.016\%$. Furthermore, we found that the mean absolute
accuracies in our astrometric calibration are $0\farcs 321$ and $0\farcs 258$
for RA and DEC, respectively. These methods would deliver even better results
for data acquired under ideal observing conditions (e.g., uniform backgrounds
and darker skies). Indeed, contrast ratios greater than $\log_{10}{(CR)}=7.73$
(i.e., 1 part in 54 million) would be easily achievable using CIDs without any
operational prerequisites.

A direct comparison between the number of cataloged sources listed on the
SIMBAD and the Gaia EDR3 databases highlights the extent to which such bright
star fields remain understudied with the Johnson-Cousins UBVRI photometric
system. Therefore, the potential next step toward practical and cost-effective
ground-based ECR imaging demonstration would be to observe different bright
star fields from a telescope with minimal internal reflections. At optimal
seeing conditions, it would be possible to detect fainter sources and achieve
even higher contrast ratios. To improve IWAs for such observations, it would
be effective to incorporate the CID imaging system with relatively simple
observing techniques, such as angular differential imaging (Marois et al.,
2006), and software-based PSF modeling techniques.

The presence of concentric circular patterns around Sirius (as seen in Figure
2 and 5) suggest the need for some development in the display tools of the
SXDR. Since the ROI is assigned during a short pre-exposure, there is a trade-
off between the pre-exposure and the atmospheric scintillation timescale. As
the SXDR continues to detect photons, the spatial positions of the brightest
pixels change within the ROI due to scintillation. To overcome this issue, a
future version of the CID technology would incorporate an option to disable
the ROI pixels by floating the gate voltages. This would indeed allow CIDs to
act as dynamic coronagraphs (Batcheldor et al., 2016).

We measured the count rate data of the detected Gaia EDR3 cataloged sources
(see Figure 2, marked in blue) to analyze the “knockdown” issue. This problem
arises due to the imager’s inability to detect photons when the pixels in the
ROI are being reset. Consequently, some flux is lost when the RAI is combined
at the end of exposure time. To characterize the knockdown coefficient, we
first calculated the contrast ratios using the cataloged magnitudes of Sirius
and the detected Gaia EDR3 cataloged sources (i.e., raw contrast ratios). We
then compared these contrast ratios to their corresponding count rate ratios.
We found that the mean knockdown coefficients for B-, V-, R-, and I-band
images are $0.91,0.92,0.90,\textrm{ and }0.92$, respectively. The observed
mean flux ratios are $28\%,29\%,24\%,\textrm{ and }33\%$ lower than their
respective true mean flux ratios.

While knockdown impedes our ability to measure the actual flux from the bright
source, this issue is beneficial under certain circumstances. For instance, it
helps increase the absolute raw contrast range possible since the light from
the bright source is unintentionally suppressed. As discussed in Batcheldor et
al. (2016) and stated above, knockdown can be easily quantified and calibrated
using the cataloged magnitudes of known sources in the bright source fields.
Therefore, for the cases where absolute photometry is required, this issue can
be resolved using the wavelength-dependent knockdown coefficients.

It is important to note that CIDs do not independently address the IWA
problem. For observations where the target PSFs are not stable or well-known,
the IWAs would be large, and detecting nearby fainter sources, such as
exoplanets, would be challenging. Since PSFs are most stable in space, CIDs
would be most beneficial on space-based telescopes. Now that the CIDs are
space-qualified to NASA’s technology readiness level 8 (Batcheldor et al.,
2020), they are excellent prospects to carry out ECR observations with future
space-based telescopes. A CID-based high-and extreme-contrast ratio astronomy
would most likely help us observe and resolve several ECR scenes, including
those with binary systems, exoplanets, circumstellar and debris disks, and
quasar host galaxies.

##  6 Conclusion

We demonstrated a direct contrast ratio of $\Delta m_{r}=18.54$,
$\log_{10}{(CR)}=7.41\pm 0.08$, or $1$ part in $26$ million using the SXDR on
the 1.0-m JKT. We implemented wavelet-based image analysis and acquired
filtered SXDR images. We applied photometric and astrometric calibration
techniques on the detected GAIA EDR3 sources and calculated contrast ratios
with respect to the cataloged magnitude of Sirius. Additionally, we detected
and resolved Sirius B at a separation of $8\farcs 401$ from Sirius. The
optical design of the JKT and high sky brightness imposed practical
limitations on the maximum achievable contrast ratios. However, direct
contrast ratios greater than $1$ part in $54$ million would be easily
achievable with CIDs at optimal observing conditions.

This research was supported in part by the American Astronomical Society’s
Small Research Grant Program and the Mount Cuba Astronomical Foundation. We
are grateful to the SARA for allocating time on the 1.0-m JKT at the Roque de
Los Muchachos Observatory, La Palma, Spain.

## References

  * 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
  * Batcheldor et al. (2016) Batcheldor, D., Foadi, R., Bahr, C., et al. 2016, Publications of the Astronomical Society of the Pacific, 128, 025001, doi: 10.1088/1538-3873/128/960/025001
  * Batcheldor et al. (2020) Batcheldor, D., Sawant, S., Jenne, J., et al. 2020, PASP, 132, 055001, doi: 10.1088/1538-3873/ab7a74
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Bhaskaran et al. (2008) Bhaskaran, S., Chapman, T., Pilon, M., & Vangorden, S. 2008, Proceedings of SPIE - The International Society for Optical Engineering, 7055, doi: 10.1117/12.795235
  * Bond et al. (2017) Bond, H. E., Schaefer, G. H., Gilliland, R. L., et al. 2017, ApJ, 840, 70, doi: 10.3847/1538-4357/aa6af8
  * Bonnet-Bidaud et al. (2000) Bonnet-Bidaud, J. M., Colas, F., & Lecacheux, J. 2000, A&A, 360, 991.  https://arxiv.org/abs/astro-ph/0010032
  * Bonnet-Bidaud & Gry (1991) Bonnet-Bidaud, J. M., & Gry, C. 1991, A&A, 252, 193
  * Bracewell & MacPhie (1979) Bracewell, R. N., & MacPhie, R. H. 1979, Icarus, 38, 136, doi: 10.1016/0019-1035(79)90093-9
  * Bradley et al. (2020) Bradley, L., Sipocz, B., Robitaille, T., et al. 2020, astropy/photutils: 1.0.0, 1.0.0, Zenodo, doi: 10.5281/zenodo.4044744
  * Eid (1995) Eid, S. I. 1995, in Electronic Imaging
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Hopkins et al. (2002) Hopkins, A. M., Miller, C. J., Connolly, A. J., et al. 2002, AJ, 123, 1086, doi: 10.1086/338316
  * Keel et al. (2016) Keel, W. C., Oswalt, T., Mack, P., et al. 2016, Publications of the Astronomical Society of the Pacific, 129, 015002, doi: 10.1088/1538-3873/129/971/015002
  * Kimble et al. (1995) Kimble, R. A., Chen, P. C., Haas, J. P., et al. 1995, in EUV, X-Ray, and Gamma-Ray Instrumentation for Astronomy VI, ed. O. H. W. Siegmund & J. V. Vallerga, Vol. 2518, International Society for Optics and Photonics (SPIE), 397 – 409, doi: 10.1117/12.218392
  * Linfield (2003) Linfield, R. P. 2003, in Interferometry in Space, ed. M. Shao, Vol. 4852, International Society for Optics and Photonics (SPIE), 443 – 450, doi: 10.1117/12.460862
  * Lowrance et al. (2005) Lowrance, P. J., Becklin, E. E., Schneider, G., et al. 2005, AJ, 130, 1845, doi: 10.1086/432839
  * Marois et al. (2006) Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556, doi: 10.1086/500401
  * Marois et al. (2008) Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348–1352, doi: 10.1126/science.1166585
  * McCreight & Goebel (1981) McCreight, C. R., & Goebel, J. H. 1981, Appl. Opt., 20, 3189, doi: 10.1364/AO.20.003189
  * McCreight et al. (1986) McCreight, C. R., McKelvey, M. E., Goebel, J. H., Anderson, G. M., & Lee, J. H. 1986, in Infrared Detectors, Sensors, and Focal Plane Arrays, ed. H. Nakamura, Vol. 0686, International Society for Optics and Photonics (SPIE), 66 – 75, doi: 10.1117/12.936526
  * Mertens & Lobanov (2015) Mertens, F., & Lobanov, A. 2015, A&A, 574, A67, doi: 10.1051/0004-6361/201424566
  * Michon, G. (1974) Michon, G. 1974, uS Patent 3,786,263
  * Miller et al. (2001) Miller, C. J., Genovese, C., Nichol, R. C., et al. 2001, AJ, 122, 3492, doi: 10.1086/324109
  * Ninkov et al. (1994) Ninkov, Z., Tang, C., Backer, B. S., Easton, R. L., & Carbone, J. M. 1994, in Astronomical Telescopes and Instrumentation
  * Oppenheimer & Hinkley (2009) Oppenheimer, B. R., & Hinkley, S. 2009, Annual Review of Astronomy and Astrophysics, 47, 253, doi: 10.1146/annurev-astro-082708-101717
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3, doi: 10.1051/0004-6361/202039587
  * Schneider et al. (2001) Schneider, G., Becklin, E. E., Smith, B. A., et al. 2001, AJ, 121, 525, doi: 10.1086/318050
  * Schneider et al. (2010) Schneider, G., Silverstone, M. D., Stobie, E., Rhee, J. H., & Hines, D. C. 2010, in Hubble after SM4. Preparing JWST, 15
  * Sparks & Ford (2002) Sparks, W. B., & Ford, H. C. 2002, ApJ, 578, 543, doi: 10.1086/342401
  * Starck & Murtagh (2002a) Starck, J.-L., & Murtagh, F. 2002a, Filtering (Berlin, Heidelberg: Springer Berlin Heidelberg), 27–58, doi: 10.1007/978-3-662-04906-8_2
  * Starck & Murtagh (2002b) —. 2002b, Introduction to Applications and Methods (Berlin, Heidelberg: Springer Berlin Heidelberg), 1–25, doi: 10.1007/978-3-662-04906-8_1
  * Starck & Murtagh (2002c) —. 2002c, Detection (Berlin, Heidelberg: Springer Berlin Heidelberg), 93–113, doi: 10.1007/978-3-662-04906-8_4
  * Starck et al. (1998) Starck, J.-L., Murtagh, F. D., & Bijaoui, A. 1998, Image Processing and Data Analysis: The Multiscale Approach (Cambridge University Press), doi: 10.1017/CBO9780511564352
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191, doi: 10.1086/131977
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9, doi: 10.1051/aas:2000332


# Spectrum-Informed Multistage Neural Network: Multiscale Function
Approximator of Machine Precision

Jakin Ng  Yongji Wang  Ching-Yao Lai

###### Abstract

Deep learning frameworks have become powerful tools for approaching scientific
problems such as turbulent flow, which has wide-ranging applications. In
practice, however, existing scientific machine learning approaches have
difficulty fitting complex, multi-scale dynamical systems to very high
precision, as required in scientific contexts. We propose using the novel
multistage neural network approach with a spectrum-informed initialization to
learn the residue from the previous stage, utilizing the spectral biases
associated with neural networks to capture high frequency features in the
residue, and successfully tackle the spectral bias of neural networks. This
approach allows the neural network to fit target functions to double floating-
point machine precision $O(10^{-16})$.

Machine Learning, ICML

##  1 Introduction

###  1.1 Precision machine learning

Typical machine learning applications, such as computer vision or natural
language processing, do not necessarily require neural networks to fit data to
extremely high precision. For instance, the training loss function may simply
be a proxy for the true metric, and moreover noise present in the training
data may cause models with very low training loss to overfit (Michaud et al.,
2023). Recently, deep learning techniques are being increasingly developed for
scientific purposes where high precision is desirable or required (Wang et
al., 2023; Michaud et al., 2023; Wang & Lai, 2024; Müller & Zeinhofer, 2023;
Jnini et al., 2024). For instance, neural networks used as interpolators for
equation discovery, which infers an exact equation or formula based on data,
the correctness of the learned equation requires very high accuracy across
multiple scales (Udrescu & Tegmark, 2020; Mojgani et al., 2024). Another
application requiring high precision is physics-informed neural networks
(PINNs), which behave as numerical solvers for partial differential equations
(PDEs), where high accuracy is an intrinsic requirement. When solving a well-
posed PDE problem, there exists a theoretical global minimum for the PINN
training where the training loss should converge to zero (Raissi et al.,
2019).

In this study, we focus on studying the regression problem as a preliminary
example to demonstrate the challenge of training neural networks to approach a
target function with high precision. The regression problem is given as: given
a dataset $\\{(\mathbf{x_{i}},u_{i}=u(\mathbf{x_{i}}))\\}$ sampled from a
continuous target function $u$, the neural network $\mathcal{N}_{\theta}$,
with parameters $\theta$, is trained to fit the data points to approximate the
function $u$. We assume that the number of data points is sufficient, and that
each of them has zero noise, which guarantees that the data set contains
sufficient authentic information of the target functions. Although the
universal approximation theorem is a theoretical guarantee that neural
networks of large enough size can approximate any function arbitrarily well
(Hornik et al., 1989), in practice, the training loss is easily trapped in
local minima and eventually plateaus after a certain number of iterations.
Many advanced techniques exist to expedite training (Sitzmann et al., 2020;
Michaud et al., 2023; Liu et al., 2020), but are unable to consistently reduce
the training error down to the required precision, which is one of the
significant challenges of using deep learning methods for scientific
objectives. In contrast, classical numerical methods can consistently reduce
error, for instance by increasing the mesh resolution.

###  1.2 Spectrum-informed initialization for regression

To resolve the precision limit of PINNs, Wang & Lai (2024) proposed the
multistage training scheme which divides neural network training into
different stages. For each stage of training, a new neural network was
introduced and optimized to learn the residues from previous stages, which
largely increases the convergence rate from linear decay to approximately
exponential decay. The combined neural networks after several stages of
training can approximate the target function up to double floating-point
machine precision $O(10^{-16})$ for one-dimensional problems. However, this
fails to hold for two- or higher-dimensional problems (Figure 1$a$). This
indicates that more advanced techniques are required to enhance the multistage
neural networks (MSNNs) to ensure high-precision approximation for higher-
dimensional problems.

The training of neural networks is known to suffer from spectral bias, which
is a phenomenon in which neural networks tend to fit low-frequency features of
the target function, and may fail to capture high-frequency information
(Rahaman et al., 2019; Xu et al., 2022). Spectral bias is particularly
problematic when using neural networks as function approximators for
multiscale problems, such as turbulence (Mojgani et al., 2024; Rybchuk et al.,
2023; Chattopadhyay & Hassanzadeh, 2023; Lai et al., 2024). Various frequency
domain approaches have been proposed to mitigate this problem, including the
scale factor approach (Cai et al., 2019; Liu et al., 2020; Jin et al., 2024;
Li et al., 2023), which were applied in the multistage neural networks (Wang &
Lai, 2024) by multiplying a large scale factor to the weight between the input
and first hidden layer. The optimal value of the scale factor was found to be
$\kappa=\pi f_{d}/\sqrt{V_{ar}}$, where $f_{d}$ is the domain frequency of the
target function and $V_{ar}$ is the variance of first layer weights. Although
this setting works effectively for fitting one-dimensional functions, it
remains limited when approximating high-frequency function in higher
dimensions. Here, to resolve the issue, we provide an advanced method that can
further alleviate spectral bias and enhance the multistage neural network
training, by initializing the neural network weights based on more
straightforward spectral information from the target function. In more
details, this spectrum-informed initialization uses the discrete Fourier
transform of the target function or the residues to inform the initialization
of the neural network, ensuring fast convergence across the spectral modes
present. In doing so, the spectrum-informed multistage neural network (SI-
MSNN) is able to approximate the target function down to $O(10^{-16})$.
Moreover, the regression provided by the neural network is accurate across the
Fourier spectrum, even for challenging multi-scale target functions.

In Section 2, we provide an overview of Fourier feature embeddings, the
discrete Fourier transform, and neural tangent kernel theory. In Section 3, we
introduce the spectrum-informed initialization of network weights, and its
usage for multi-stage neural networks (MSNNs). In Section 4, we provide
experimental results demonstrating the advantage of SI-MSNNs over the original
MSNNs on error reduction and its ability to reach machine precision for
complex regression problems. Lastly, in Section 5, we provide further
discussions and future work on the spectrum-informed initialization for
multistage neural networks.

Figure 1: ($a$) Target function $u_{g}(x,y)$ for a 2-D regression problem.
($b$) Errors of neural networks with the target function $u_{g}$ after
different stages of training using the original multi-stage neural networks
(MSNNs) where the scaling factor was used to mitigate the spectral biases of
the network training. The 3rd-stage residue reaches $O(10^{-8})$. ($c$) Errors
after different stages of training using spectrum-informed initialization for
both weights and biases of the network, which reaches $O(10^{-13})$ by the 3rd
stage of training.

##  2 Preliminaries

###  2.1 Discrete Fourier transform

The Fourier transform decomposes a function into the various frequencies and
corresponding amplitudes that are present, writing the function in terms of a
complex exponential or sinusoidal/cosinusoidal basis (Sneddon, 1995). The
discrete Fourier transform (DFT) of a function with values provided at $N$
equally spaced points can be computed using the fast Fourier transform (FFT)
algorithm in $O(N\log N)$. The one-dimensional DFT of a function $u$,
specified at $N$ points $x_{n}$, is given by

$\tilde{u}(k)=\frac{1}{N}\sum_{n=0}^{N-1}u(x_{n})\exp(-ikx_{n}),$ (1)

for which the inverse Fourier transform is

$\tilde{u}(k)=\sum_{n=0}^{N-1}u(x_{n})\exp(ikx_{n}),$ (2)

where $k$ is the wavenumber. Analogously, the two-dimensional DFT of a
function specified on a grid with $N_{x}$ values $x_{n}$ and $N_{y}$ values
$y_{m}$ is given by

$\tilde{u}(\mathbf{k})=\sum_{n=0}^{N_{x}}\sum_{m=0}^{N_{y}}\exp(-i(k_{x}x_{n}+k_{y}y_{m})),$
(3)

where $\mathbf{k}=(k_{x},k_{y})$ is the wavenumber vector.

###  2.2 Fourier feature networks

Spectral bias, which is also referred to as the frequency principle, is a
well-known phenomenon of standard multi-layer perceptrons (MLPs), which tend
to exhibit a bias towards fitting low-frequency information in target
functions, and learn high-frequency information more slowly, or fail to
converge on high frequency modes (Rahaman et al., 2019; Xu et al., 2022). To
allow MLPs to learn high-frequencies, Fourier feature networks first pass
input points $x$ through a Fourier feature mapping $\gamma(x)$, before passing
the result through the MLP (Rahimi & Recht, 2007; Tancik et al., 2020). The
Fourier feature mapping $\gamma(x)$ can be alternatively considered as the
first layer of the MLP, as a composition of a nonlinear, periodic activation
function and a linear function with added bias weights.

####  2.2.1 Random Fourier features

Random Fourier features are a non-trainable input mapping

$\gamma(x)=\cos(2\pi Bx+b)$ (4)

where $B$ and $b$ are fixed random variables, and the cosinusoidal function is
applied component-wise (Rahimi & Recht, 2007).

####  2.2.2 Fourier feature mapping

Fourier feature mapping extends this technique, using the positional encoding

$\gamma(x)=[A\cos(2\pi Bx),A\sin(2\pi Bx)]^{T},$ (5)

where $B\in\mathbb{R}^{m\times d}$ are the weights for the feature mapping,
which can be either fixed or trainable. Typically, these weights are
initialized as $B\sim\mathcal{N}(0,\sigma^{2})$, where $\sigma$ is a
hyperparameter that determines the span of the initialized weights and can be
set to 1 as a default, or determined based on the problem (Tancik et al.,
2020).

####  2.2.3 Cosinusoidal activation

Since cosine and sine differ only by a phase, which can be recovered using the
trainable phase parameter $b$, it is comparable to use the cosinusoidal
mapping

$\gamma(x)=[A\cos(2\pi Bx+b)]^{T},$ (6)

where the weights $B$ are initialized as in the previous section, and the
biases $b$ can be initialized to zero or following a normal distribution.
Equivalently, the cosinusoidal feature mapping can be seen as a modified first
layer of an MLP, where the activation function is $\sigma(x)=A\cos(2\pi x)$
rather than the typical choices $\sigma(x)=\tanh(x)$ or
$\sigma(x)=\operatorname{ReLU}(x).$

###  2.3 Neural tangent kernel

The neural tangent kernel (NTK) theory describes the evolution of neural
networks by gradient descent during training. In the limit of an infinite-
width neural network, the NTK can be used to approximate the result of
training a neural network as the learning rate approaches zero. In fact, the
neural network converges in the eigenvectors of the NTK at an exponential rate
with respect to the corresponding eigenvalue (Jacot et al., 2018; Tancik et
al., 2020). For a conventional MLP, the spectral bias can be viewed as an
eigenvector bias, as the eigenvalues decay rapidly (Wang et al., 2021). When
using a Fourier feature mapping, the principal eigenvectors of the NTK
correspond to the embedded frequencies from $B$, and thus embedding the
frequencies of the target function $u$ into the neural network using Fourier
feature mapping allows rapid convergence in those directions (Tancik et al.,
2020).

##  3 Methods

###  3.1 Problem setup

We consider the supervised learning regression problem, where $N_{d}$ data
points $\\{(\mathbf{x_{i}},u(\mathbf{x_{i}}))\\}$ are provided, and the task
is to fit a neural network to the target function $u$. In the precision
machine learning setting, we assume that the data is exact, in order to test
the capacity of neural networks as universal function approximators.

A neural network with $L$ hidden layers, given a feature embedding
$\gamma(x)$, is a mapping $\mathcal{N}:\mathbb{R}^{d}\to\mathbb{R}$ which acts
as
$\mathcal{N}(x;\theta)=f_{L}\circ\sigma_{L}\circ\cdots\cdots\circ\sigma_{1}\circ
f_{1}\circ\gamma(x)$, where $f_{k}(x)=W_{k}x+b_{k}$ is a linear function and
$\sigma_{k}(x)$ is a nonlinear activation function, such as $\tanh$ or
$\operatorname{reLU}$.

###  3.2 Multi-stage neural networks (MSNN)

As mentioned above, multi-stage neural networks (Wang & Lai, 2024) are a novel
approach dividing the training into stages, maintaining a high rate of
convergence and allowing the regression for one-dimensional problems to reach
machine precision (Wang & Lai, 2024).

Given a target function $u$, a typical neural network $u_{0}(x)$ with Xavier
weight initialization is trained to approximate $u/\epsilon_{0}$, with the
input coordinates $\mathbf{x}$ normalized to within $[-1,1]$. The normalizing
factor $\epsilon_{0}$ is taken to be the root mean square value
$\epsilon_{0}=\sqrt{\frac{1}{N_{d}}\sum_{i=0}^{N_{d}-1}u_{i}^{2}}$, where
$u/\epsilon_{0}$ has order of magnitude $O(1)$ matching the output of the
network at initialization. The first stage residue is
$e_{1}(x)=u(x)-u_{0}(x)$, and is typically a high-frequency function.

Then, the second-stage neural network is trained to approximate the normalized
first stage residue $e_{1}/\epsilon_{1}$ as a target function, where
$\epsilon_{1}$ is the root mean square value of the error $e_{1}$ on the
provided data points. To effectively fit the high-frequency first stage
residue given the dominant frequency $f_{d}$ and the variance of weights in
the first layer $V_{ar}$, the activation function of the first layer is taken
to be the periodic, high-frequency $\sigma(x)=\cos(\kappa x)$, where its
inputs are first multiplied by the optimal scale factor $\kappa=\pi
f_{d}/\sqrt{V_{ar}}$ (Wang & Lai, 2024). The re-scaling of the input to enable
learning of high-frequency functions has also been implemented as adaptive
activation functions (Jagtap et al., 2020).

Figure 2: (a) Target function $\psi(x,y)$, a single time snapshot of the
numerical solution of the stream function for the 2D incompressible Navier-
Stokes equations (11) and (12) with Re = 2000. (b–e) The residues after each
of the four stages of training, which are the target functions for the next
stages. After four stages, the residue has approached machine precision
$O(10^{-16})$. (g) Result of the Spectrum-Informed Multistage Neural Network
(SI-MSNN) with four stages of training. (h–j) The spectral domain of each of
the residues.

In general, the $(n+1)$th stage neural network $u_{n+1}$ is the previous stage
residue $e_{n}(x)=u(x)-\sum_{i=1}^{n}\varepsilon_{n}u_{n}(x)$, normalized by
its root mean square value $\epsilon_{n}$. The final result with $s$ stages is
$\sum_{n=1}^{s}\varepsilon_{n}u_{n}(x)\approx u(x).$ By tuning the neural
network in each stage to capture the desired magnitudes and frequencies of the
residue, the multistage neural network approach is able to reach machine
precision $O(10^{-16})$ for the regression problem in one dimension.

###  3.3 $L^{p}$ Norm

The typical loss function for training and evaluation in regression problems
is the mean squared error (MSE), which is the squared $L^{2}$ norm, defined as
$\|x\|_{L_{2}^{2}}=\sum_{n=1}^{N_{d}}x_{n}^{2}.$ A more general loss function
uses the $L^{p}$ norm, where

$\|x\|_{L_{p}}=\left(\sum_{n=1}^{N_{d}}|x_{n}|^{p}\right)^{1/p},$ (7)

which has a comparable order of magnitude to the root mean squared error
(RMSE). As $p$ increases, larger deviations are more heavily penalized. For
our purposes, we use the $L^{p}$ norm with $p=10$, which prevents large spikes
in the residue so the next stage of the multistage neural network can learn
the residue more readily.

###  3.4 Spectrum-informed initialization

In two or more dimensions, the multistage neural network is not able to reach
machine precision, and fails to capture the high frequencies present in the
later stage residues. We propose a spectrum-informed initialization to replace
the cosinusoidal mapping with a scale factor $\kappa$ used in Wang & Lai
(2024). Instead, here we tailor the neural network to the specific frequencies
present in the dataset. Consider a dataset with input-output pairs
$(\mathbf{x}_{i},u_{i}=u(\mathbf{x}_{i}))$ for $\mathbf{x}_{i}\in[-1,1]^{2}$
on a grid with $N$ points in both dimensions, where the input domain is
assumed to already be transformed to $[-1,1]$, and $u_{i}\in\mathbb{R}$. The
spectrum-informed initialization provides a targeted initialization of $B$ and
$b$ for the input mapping layer $\gamma(x;B,b)=[A\cos(Bx+b)]^{T}$ of a neural
network predicting the quantity $u(x)$ given the input coordinates $x$, where
$B$ and $b$ are tunable parameters and $A$ is a fixed hyperparameter. The
remainder of the neural network layers can be initialized with the usual
Xavier initialization scheme, which prevents gradient vanishing or explosion
(Glorot & Bengio, 2010).

Algorithm 1 Spectrum-informed initialization

Compute the discrete Fourier transform $\tilde{u}$ of $u$, where
$\tilde{u}(\mathbf{k})=a(\mathbf{k})e^{i\theta(\mathbf{k})}$ in polar form.
Let $\alpha(\mathbf{k})=\frac{2}{N^{2}}a(\mathbf{k})$ if $k_{y}>0$, and
$\alpha(\mathbf{k})=\frac{1}{N^{2}}a(\mathbf{k})$ if $k_{y}=0$ be a normalized
version of the magnitude.

Let $\mathbf{k}^{(j)}$ be the Fourier mode corresponding to the $j$th largest
magnitude $|\tilde{u}(\mathbf{k}^{(j)})|$, considering only the modes such
that $k_{y}\geq 0.$ Take $\alpha^{(j)}=\alpha(\mathbf{k}^{(j)})$ and
$\theta^{(j)}=\theta(\mathbf{k}^{(j)})$.

Given a desired first layer width $n_{f}$, let
$B=[\mathbf{k}^{(1)},\cdots,\mathbf{k}^{(n_{f})}]^{T}$,
$b=[\theta^{(1)},\cdots,\theta^{(n_{f})}]^{T}$ and
$A=\left[\alpha^{(1)},\cdots,\alpha^{(n_{f})}\right]^{T}$

Initialize the first layer of the neural network, which has width $n_{f}$,
with the Fourier feature mapping $\gamma(x;B,b)=[A\cos(Bx+b)]^{T}$. Initialize
the other parameters based on the Xavier scheme.

The two-dimensional DFT of $u$ satisfies

$\displaystyle u(x,y)$
$\displaystyle=\frac{1}{N^{2}}\sum_{k_{x},k_{y}}\tilde{u}(k_{x},k_{y})\exp(i(k_{x}x+k_{y}y))$
(8)
$\displaystyle=\sum_{k_{x},k_{y}}\frac{1}{N^{2}}a(\mathbf{k})\exp(i(\mathbf{k}\cdot\mathbf{x}+\theta(\mathbf{k})))$
(9) $\displaystyle=\sum_{k_{x},k_{y}\geq
0}\alpha(\mathbf{k})\cos(\mathbf{k}\cdot\mathbf{x}+\theta(\mathbf{k})),$ (10)

where $\alpha({\bf k})$ are as defined in Algorithm 1. The last equality comes
from the conjugate symmetry
$\tilde{u}(\mathbf{k})=\overline{\tilde{u}(-\mathbf{k})}$ of $\tilde{u}$,
since $u$ is a real-valued function. Using the conjugate symmetry avoids
embedding redundant information.

Then, given the $n_{f}$ largest modes, the Fourier feature mapping, which can
be considered the first layer of the neural network, is
$\gamma(x)=A\cos(Bx+b)=\sum_{i\leq
n_{f}}\alpha^{(i)}\cos(\mathbf{k}^{(i)}\cdot\mathbf{x}+\theta^{(i)})$, which
contains precisely the information from the top $n_{f}$ modes of the Fourier
transform of $u$. If $n_{f}=N^{2}/2$ and the neural network contains only one
layer, which is the Fourier feature embedding, the single-layer neural network
will exactly match the DFT representation of $u$, as in Equation 10.

Based on the neural tangent kernel (NTK) theory (Jacot et al., 2018),
embedding these primary Fourier modes $\mathbf{k}$ allows the neural network
to effectively learn the target function across the frequency spectrum.

##  4 Experiments

In this section, we experimentally demonstrate the ability of the spectrum-
informed initialization to learn the target frequencies, and to train neural
networks as approximators down to the limits of numerical precision.

###  4.1 Comparison with scale factor approach

In Figure 1, we compare the results of a three-stage neural net with the
original scale factor approach and the spectrum-informed initialization,
fitting a simple target function $u_{g}(x)=[\sin(2x+1)-0.5x](1-y^{2})$. For
the scale factor approach, the Fourier feature embedding is of the form
$\gamma(x)=\cos(\kappa w^{(0)}x+b),$ where $w^{(0)}$ are the first-layer
weights following the Gaussian distribution $\mathcal{N}(0,V_{ar})$ with the
variance $V_{ar}$ set by the Xavier initialization method. The scale factor
$\kappa=\pi f_{d}/\sqrt{V_{ar}}$ is determined based on the dominant frequency
$f_{d}$ of the target function. In contrast, the spectrum-informed
initialization employs the Fourier feature embedding $\gamma(x)=A\cos(Bx+b)$
with $A,B,$ and $b$ being directly initialized by the spectral information of
the target function, as described in Algorithm 1.

Experimentally, the scale factor approach yields a 3rd-stage residue of
$O(10^{-8})$ (Figure 1$b$), whereas the spectrum-informed initialization
allows the MSNN to fit the data up to $O(10^{-13})$ (Figure 1$c$), which is
five orders of magnitudes more accurate. We note that the number of iterations
used in each stage of the two experiments are the same.

Figure 3: Left: A comparison of the power spectrum of the 2D Navier-Stokes
example in Fig. 2 given by a single-stage SI-MSNN and a four-stage SI-MSNN,
each with three hidden layers of width 30. The single-stage SI-MSNN has first
layer width $n_{f}=10000$. The first stage of the four-stage SI-MSNN has first
layer width $n_{f}=1359$, based on the number of primary Fourier modes
present, and the remaining stages all have first layer width $n_{f}=10000.$
Right: The loss convergence of a single-stage SI-MSNN compared to a four-stage
SI-MSNN, where training terminates when machine precision is reached.

###  4.2 2D incompressible Navier-Stokes

As an example problem, we consider fitting a snapshot of two-dimensional
homogeneous isotropic decaying turbulence (2D-DHIT), which is a broadly multi-
scale problem. Random Fourier features have been shown to be effective for
similar multi-scale turbulence problems (Mojgani et al., 2024). The Fourier
coefficients of the target function decay exponentially with the wavenumber,
where the high-frequency, small-scale information plays a large role in the
behavior of the system. Taking the curl of the 2D incompressible Navier-Stokes
equations yields the vorticity equation, which can be written in terms of the
stream function $\psi$ and the vorticity $\omega$:

$\frac{\partial\omega}{\partial t}+\frac{\partial\psi}{\partial
y}\frac{\partial\omega}{\partial x}-\frac{\partial\psi}{\partial
x}\frac{\partial\omega}{\partial y}=\frac{1}{\text{Re}}\nabla^{2}\omega$ (11)

and

$\nabla^{2}\psi=-\omega,$ (12)

which can model a variety of geophysical flows (Subel et al., 2023). Using a
direct numerical simulation (DNS), on a $[0,2\pi)^{2}$ grid with resolution
$N=512$ and Reynolds number $Re=2000$, the 2D-DHIT equations are solved using
a pseudo-spectral solver with third-order Runge-Kutta time stepping, with the
timestep determined by the CFL condition. The initial condition is taken to be
a random vorticity field with a broad banded energy spectrum (McWilliams,
1984).

We perform regression with the target function as a snapshot of the stream
function $\psi(\mathbf{x},t=1.25)$. The coordinate inputs to the neural
network are re-scaled to the domain $[-1,1]$, and the spectrum-informed
initialization is based on the re-scaled coordinates. The results are
transformed back to the original coordinates from $[0,2\pi)$. In terms of
experimental setup, we employ a fully-connected neural network with the first
layer as the spectrum-informed Fourier feature embedding with cosine as the
activation function, as well as three hidden layers with 30 units in each
layer using hyperbolic tangent activation functions. The first layer contains
$n_{f}=1359$ units for the first-stage neural network, and $n_{f}=10000$ units
for the second, third, and fourth-stage neural networks, based on the number
of primary Fourier modes in the target function of each stage. The spectrum-
informed initialization is used for the first Fourier feature layer, and the
Xavier scheme is used to initialize the other parameters, as described in
Algorithm 1. The neural networks are trained with full batch gradient descent
using Adam optimizer with initial learning rate of $10^{-3}$ and incorporating
learning rate annealing. The models are trained for 150000 iterations, with
the fourth stage terminating when machine precision is reached.

In Figure 2, the SI-MSNN method is successfully able to reduce the error down
to the machine precision of a 64-bit double float, $O(10^{-16})$, after four
stages of training. This result experimentally validates the ability of SI-
MSNNs to approximate a target function with arbitrarily high accuracy.

In Figure 3, the multistage setup is shown to be necessary. With a single
stage of training, the loss convergence plateaus to a linear decay rate, and
the neural network is not able to correct the residues relative to the target
function. However, by incorporating multiple stages, with a neural network
initialized based on the magnitude and spectrum of the residue, an exponential
convergence rate can be maintained, allowing efficient learning.

Moreover, the spectrum-informed initialization allows the multistage neural
network to accurately represent the function across the spectral domain. Here,
the two-dimensional Fourier transform of the stream function $\psi$ satisfies
$\tilde{\psi}(k_{x},k_{y})=\sum_{n=0}^{N_{x}}\sum_{m=0}^{N_{y}}\psi(x_{n},y_{m})\exp(-i(k_{x}x_{n}+k_{y}y_{m}))$,
where $\mathbf{k}=(k_{x},k_{y})$ is the wavenumber vector. Given the Fourier
transform $\tilde{\psi}$, the angle-averaged power spectrum $E(k)$ is defined
as

$E(k)=\sum_{k-\Delta k\leq|\mathbf{k}|\leq k+\Delta
k}|\tilde{\psi}(\mathbf{k})|^{2},$ (13)

for a given spectral band $\Delta k$ (Kag et al., 2022). The spectral band is
defined so that the wavenumber $k$ falls into 200 bins.

In Figure 3, as the number of stages increases, the SI-MSNN is able to match
the power spectrum $E(k)$ of the target function $\psi(x,y)$ up to increasing
wavenumbers. Successfully fitting both the target function and its power
spectrum is important for multi-scale scientific problems, and the capacity to
do so is shown with this 2D Navier-Stokes example.

##  5 Discussions

We introduce a novel spectrum-informed initialization, allowing efficient
training of neural networks for solving the regression problem to machine
precision. By utilizing the spectral biases of neural networks, a spectrum-
informed Fourier embedding of the input allows the neural network to learn in
the spectral domain and converge rapidly. The spectrum-informed initialization
for multistage neural networks (SI-MSNN) allows the neural network to fit
target functions down to machine precision, even for multi-scale target
functions, as validated experimentally. The results using the SI-MSNN
demonstrate that this approach can achieve residues many orders of magnitude
smaller than state-of-the-art approaches, down to machine precision.

In the future, we plan on applying the spectrum-informed initialization for
scientific machine learning problems requiring precision as high as possible.
In particular, we propose using the spectrum-informed initialization for
multi-stage physics-informed neural networks as partial differential equation
solvers (Wang & Lai, 2024; Raissi et al., 2019), where the boundary and
initial conditions are the error-free data provided, along with a set of
governing equations. Our work provides a promising approach towards precision
machine learning.

## References

  * Cai et al. (2019) Cai, W., Li, X., and Liu, L.  Phasednn-a parallel phase shift deep neural network for adaptive wideband learning.  _arXiv preprint arXiv:1905.01389_ , 2019.
  * Chattopadhyay & Hassanzadeh (2023) Chattopadhyay, A. and Hassanzadeh, P.  Long-term instabilities of deep learning-based digital twins of the climate system: The cause and a solution.  _arXiv preprint arXiv:2304.07029_ , 2023.
  * Glorot & Bengio (2010) Glorot, X. and Bengio, Y.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pp. 249–256. JMLR Workshop and Conference Proceedings, 2010.
  * Hornik et al. (1989) Hornik, K., Stinchcombe, M., and White, H.  Multilayer feedforward networks are universal approximators.  _Neural networks_ , 2(5):359–366, 1989.
  * Jacot et al. (2018) Jacot, A., Gabriel, F., and Hongler, C.  Neural tangent kernel: Convergence and generalization in neural networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Jagtap et al. (2020) Jagtap, A. D., Kawaguchi, K., and Karniadakis, G. E.  Adaptive activation functions accelerate convergence in deep and physics-informed neural networks.  _Journal of Computational Physics_ , 404:109136, 2020.
  * Jin et al. (2024) Jin, G., Wong, J. C., Gupta, A., Li, S., and Ong, Y.-S.  Fourier warm start for physics-informed neural networks.  _Engineering Applications of Artificial Intelligence_ , 132:107887, 2024.
  * Jnini et al. (2024) Jnini, A., Vella, F., and Zeinhofer, M.  Gauss-newton natural gradient descent for physics-informed computational fluid dynamics.  _arXiv preprint arXiv:2402.10680_ , 2024.
  * Kag et al. (2022) Kag, V., Seshasayanan, K., and Gopinath, V.  Physics-informed data based neural networks for two-dimensional turbulence.  _Physics of Fluids_ , 34(5), 2022.
  * Lai et al. (2024) Lai, C.-Y., Hassanzadeh, P., Sheshadri, A., Sonnewald, M., Ferrari, R., and Balaji, V.  Machine learning for climate physics and simulations.  _arXiv preprint arXiv:2404.13227_ , 2024.
  * Li et al. (2023) Li, S., Xia, Y., Liu, Y., and Liao, Q.  A deep domain decomposition method based on fourier features.  _Journal of Computational and Applied Mathematics_ , 423:114963, 2023.
  * Liu et al. (2020) Liu, Z., Cai, W., and Xu, Z.-Q. J.  Multi-scale deep neural network (mscalednn) for solving poisson-boltzmann equation in complex domains.  _arXiv preprint arXiv:2007.11207_ , 2020.
  * McWilliams (1984) McWilliams, J. C.  The emergence of isolated coherent vortices in turbulent flow.  _Journal of Fluid Mechanics_ , 146:21–43, 1984.
  * Michaud et al. (2023) Michaud, E. J., Liu, Z., and Tegmark, M.  Precision machine learning.  _Entropy_ , 25(1):175, 2023.
  * Mojgani et al. (2024) Mojgani, R., Chattopadhyay, A., and Hassanzadeh, P.  Interpretable structural model error discovery from sparse assimilation increments using spectral bias-reduced neural networks: A quasi-geostrophic turbulence test case.  _Journal of Advances in Modeling Earth Systems_ , 16(3):e2023MS004033, 2024.
  * Müller & Zeinhofer (2023) Müller, J. and Zeinhofer, M.  Achieving high accuracy with pinns via energy natural gradients.  _arXiv preprint arXiv:2302.13163_ , 2023.
  * Rahaman et al. (2019) Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F., Bengio, Y., and Courville, A.  On the spectral bias of neural networks.  In _International conference on machine learning_ , pp. 5301–5310. PMLR, 2019.
  * Rahimi & Recht (2007) Rahimi, A. and Recht, B.  Random features for large-scale kernel machines.  _Advances in neural information processing systems_ , 20, 2007.
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., and Karniadakis, G. E.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational physics_ , 378:686–707, 2019.
  * Rybchuk et al. (2023) Rybchuk, A., Hassanaly, M., Hamilton, N., Doubrawa, P., Fulton, M. J., and Martínez-Tossas, L. A.  Ensemble flow reconstruction in the atmospheric boundary layer from spatially limited measurements through latent diffusion models.  _Physics of Fluids_ , 35(12), 2023.
  * Sitzmann et al. (2020) Sitzmann, V., Martel, J., Bergman, A., Lindell, D., and Wetzstein, G.  Implicit neural representations with periodic activation functions.  _Advances in neural information processing systems_ , 33:7462–7473, 2020.
  * Sneddon (1995) Sneddon, I. N.  _Fourier transforms_.  Courier Corporation, 1995.
  * Subel et al. (2023) Subel, A., Guan, Y., Chattopadhyay, A., and Hassanzadeh, P.  Explaining the physics of transfer learning in data-driven turbulence modeling.  _PNAS nexus_ , 2(3):pgad015, 2023.
  * Tancik et al. (2020) Tancik, M., Srinivasan, P., Mildenhall, B., Fridovich-Keil, S., Raghavan, N., Singhal, U., Ramamoorthi, R., Barron, J., and Ng, R.  Fourier features let networks learn high frequency functions in low dimensional domains.  _Advances in neural information processing systems_ , 33:7537–7547, 2020.
  * Udrescu & Tegmark (2020) Udrescu, S.-M. and Tegmark, M.  Ai feynman: A physics-inspired method for symbolic regression.  _Science Advances_ , 6(16):eaay2631, 2020.
  * Wang et al. (2021) Wang, S., Wang, H., and Perdikaris, P.  On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics-informed neural networks.  _Computer Methods in Applied Mechanics and Engineering_ , 384:113938, 2021.
  * Wang & Lai (2024) Wang, Y. and Lai, C.-Y.  Multi-stage neural networks: Function approximator of machine precision.  _Journal of Computational Physics_ , pp. 112865, 2024.
  * Wang et al. (2023) Wang, Y., Lai, C.-Y., Gómez-Serrano, J., and Buckmaster, T.  Asymptotic self-similar blow-up profile for three-dimensional axisymmetric euler equations using neural networks.  _Physical Review Letters_ , 130(24):244002, 2023.
  * Xu et al. (2022) Xu, Z.-Q. J., Zhang, Y., and Luo, T.  Overview frequency principle/spectral bias in deep learning.  _arXiv preprint arXiv:2201.07395_ , 2022.