11institutetext: Institut für Experimentelle und Angewandte Physik,
Universität Kiel, Kiel, Germany 11email<EMAIL_ADDRESS>22institutetext: Universidad de Alcalá, Space Research Group, 28805 Alcalá de
Henares, Spain 33institutetext: National Observatory of Athens/IAASARS,
Athens, Greece 44institutetext: SMD/Heliophysics Division, NASA HQ,
Washington, DC, USA, 55institutetext: Max-Planck-Institut für
Sonnensystemforschung, Göttingen, Germany.

# The Electron Proton Helium INstrument as an Example for a Space Weather
Radiation Instrument

Patrick Kühl 11 Bernd Heber 11 Raúl Gómez-Herrero 22 Olga Malandraki 33 Arik
Posner 44 Holger Sierks 55

###### Abstract

The near-Earth energetic particle environment has been monitored since the
1970’s. With the increasing importance of quantifying the radiation risk for,
e.g. for the human exploration of the Moon and Mars, it is essential to
continue and further improve these measurements. The Electron Proton Helium
INstrument (EPHIN) on-board SOHO continually provides these data sets to the
solar science and space weather communities since 1995. Here, we introduce the
numerous data products developed over the years and present space weather
related applications. Important design features that have led to EPHINs
success as well as lessons learned and possible improvements to the instrument
are also discussed with respect to the next generation of particle detectors.

###### keywords:

energetic particle detector – space weather instrumentation – cosmic rays,
solar energetic particle events

##  1 Introduction

Solar energetic particles, from suprathermal (few keV) up to relativistic (few
GeV for protons and ions) energies are an important contributor to the space
environment characterization. They are emitted from the Sun in association
with solar flares and Coronal Mass Ejection (CME)-driven shock waves. SEP
events constitute a serious radiation hazard and pose a threat to modern
technology which relies on spacecraft as well as to humans in space. In
addition they are of concern for avionics and commercial aviation. Thus,
mitigation procedures must be developed. Novel SEP event forecasting tools
have been developed within the HESPERIA H2020 EU project and are strongly
relied upon to mitigate against SEP events.

Figure 1: Sketch (left) and photograph (right) of the EPHIN instrument.

These forecasting tools as well as the scientific studies of the very events
they are designed for to forecast, naturally share some limitations such as
the availability and quality of the underlying data. Arguably one of the most
import data sources for space weather applications is the NASA/ESA SOlar and
Heliospheric Observatory (SOHO) that was launched in 1995 and is since 1996
orbiting the Lagrangian point L1. The scientific payload of the spacecraft
consists of several remote and in-situ instruments including the Electron
Proton Helium Instrument (EPHIN), a particle telescope with a field of view of
about 83∘ and a geometry factor of 5.1 cm2 sr that measures electrons with
energies between 0.25 and 10.4 MeV as well as protons and Helium in the energy
range of 4.3 up to above 53 MeV/nucleon (Müller-Mellin et al., 1995). A
schematic view of the EPHIN instrument is shown in Figure 1. The instrument
consists of a stack of six Solid State Detectors (labeled A-F in Fig. 1) that
are surrounded with a scintillation detector (G), that acts as an anti-
coincidence. Stopping ions are identified by applying the
$\frac{dE}{dx}-E$-method. Since the telemetry data rate is limited, energy
losses in SSDs A-E can only be transmitted for a statistical ensemble of all
particles counted by EPHIN. This constitutes the so-called Pulse Height
Analysis (PHA) data. A crude onboard separation between electrons, protons and
helium has been introduced for different penetration depths, i.e. different
coincidences, by comparing the energy-loss in SSD A to different thresholds.
Futhermore, SSDs A and B are segmented (Fig. 1) with the outer segments being
switched-off automatically by the instrument during periods of high fluxes
(i.e. if the count rate of the central segment of detector A is above 15000
counts per second) in order to reduce the geometric factor by a factor of 24
resulting in fewer electronic pile-up and dead time issues (for details see
Müller-Mellin et al., 1995).
The instrument has been proven to be of significant value not only for the
space science but also for the space weather community especially due to its
consistent measurements over two and a half decades as well as its high alive
time during the mission. Certain long-term measurements such as those of MeV
electrons with a low background are even unique to EPHIN. Apart from the
communication loss to SOHO in 1998, EPHIN has a data coverage of over 90%.
Figure 2 clearly shows that most of the data gaps correspond with roll
maneuvers of the SOHO spacecraft (indicated by blue lines). The reason for the
SOHO roll maneuvers are mechanical issues with SOHO’s high gain antenna
starting in 2003 (Policella et al., 2009). In order to reduce the limitations
in telemetry that come with an antenna unable to actively point to Earth, it
was decided to move the antenna one last time to a ”sweet spot” that could
guarantee a stable connection to Earth during half of SOHOs orbit around L1
and rotate the entire SOHO spacecraft twice per orbit (i.e. every three
months). In-between both half-orbits (around the time of the roll maneuver),
however, a short period during which communication is limited remains. These
keyhole periods are the reason for the majority of EPHIN’s data gaps. Note
that the pointing of EPHIN (nominal: 45 deg west from sun-s/c line, i.e. along
the nominal Parker spiral) is also affected by the roll maneuvers (pointing
becomes perpendicular to the nominal spiral). The dates of the roll maneuvers
as well as the attitude of spacecraft can be found here:
https://soho.nascom.nasa.gov/data/ancillary/attitude/roll/nominal_roll_attitude.dat.
The keyhole periods are also listed here:
https://sohowww.nascom.nasa.gov/soc/keyholes.txt
Despite being active and measuring particles for almost 25 years - the primary
mission duration was defined as two years (Domingo et al., 1995) -, forced
inactivity and possible harsh temperature environments for three months during
the communication loss of SOHO in 1998 (Vandenbussche, 2006) as well as being
exposed to various large SEP events (Kühl and Heber, 2019, and references
therein), EPHIN remains in a remarkably good condition and is still providing
data for and fulfilling most of its primary scientific objectives. The
negative impact of the only two technical issues, increased noise level 1) in
detector E in 1997 and in 2) in detector D in 2017, was limited by careful
software management and data analysis. After the noise increases in the
detectors were monitored, EPHIN was commanded to failure modes ”FME” and ”FMD”
on February, 19th, 1997 and October, 4th, 2017, respectively. The failure
modes remove the defective detectors from the coincidence logic which results
in lower number of energy channels per particle type (four in nominal mode,
three and two for FME and FMD, respectively). It was possible, however, to
retrieve the nominal energy channels based on the energy losses in the SSDs A,
B and C using PHA data (see section 2.3).
Section 2 provides an overview of all different data products of EPHIN
including advantages and disadvantages as well as applications of the
different data products. Section 3 describes applications of the EPHIN data
with respect to space weather such as investigations and forecasts of SEP
events, in section 4 limitations and possible improvements for EPHIN are
discussed based on the lessons learned. A summary concludes the importance of
measurements from EPHIN for the space weather community and proposes design
features for future instruments.

Figure 2: Alive time of EPHIN over six years. The blue lines indicate the date
of a SOHO roll maneuver.

##  2 EPHIN Data Products

The raw (level 0) EPHIN data telemetered to ground is processed and formatted
routinely at the Christian-Albrechts University of Kiel. Several data levels
and file formats have been implemented and released publicly over the years.
All data products are provided as merged daily files with a cadence of 1
minute except for the histogram data (see section 2.1.2), which has a time
resolution of eight minutes. The following section introduces the most
important data products, flaws and advantages, download sources for the data
as well as their usage for space weather applications. While the majority of
these data sets are available only with a couple of days of delay after
acquisition onboard, some limited data are also available in near real-time
(see section 3.2). Documentations of the different data products and their
file format can be found here: http://ulysses.physik.uni-kiel.de/costep/doc/

###  2.1 Level 1

The Level 1 data products consist of 1) housekeeping data (not covered here),
2) single & coincidence counter as well as 3) histograms, and 4) PHA data. The
SCI files that include the counters and histograms as well as the PHA files
can be downloaded from http://ulysses.physik.uni-kiel.de/costep/level1/

####  2.1.1 Single Counter Rates

Independent of any coincidence logic, EPHIN tracks the count rate of each
individual detector (including the anticoincidence G). Whether or not a
particle is counted in a detector depends on its energy deposition in the
detector in comparison to the detectors lower threshold (energy thresholds are
given in table 3, Müller-Mellin et al., 1995). While originally included in
the data stream in order to monitor the noise level of the detectors, it has
been shown that these single counter rates are also of scientific value due to
their high statistics (an idea first proposed by Richardson et al., 1996).
Heber et al. (2015) has shown that EPHIN is capable of measuring count rate
variations with amplitudes in the below 1% regime on an hourly time resolution
and hence, observe Forbush Decreases associated to ICMEs in general and
stealth CMEs in particular (see for example Howard and Harrison, 2013).
Utilizing the anti-coincidence detector G the measurement capabilities allow
to investigate FDs using accumulation times as low as 10 minutes as shown in
Figure 3.

Figure 3: Variation of 10 minute averaged count rates of the single detector
rate G (red curve) and the integral channel (black curve) during the passage
of a Stealth CME (Heber et al., 2015). For details see text.

The figure shows a time period discussed in detail in Heber et al. (2015) when
a Forbush Decrease (FD) caused by Stealth CME was observed. It compares the
count rate variations of the single detector G (red curve) to that of the
integral channel (i.e. coincidence of SSDs A to F, black curve). Both data
sets have been averaged over 10 minutes and are normalized to the count rate
measured during the time period from 19:00 to 22:30 UT on July 20, 2009. While
the count rate for the single detector is about 27,000 counts/minute the
integral channel counts about 240 counts/minute leading to a 10 times higher
statistical accuracy of the single detector compared to the integral count
rate. The corresponding statistical uncertainty of the measurements shown in
Figure 3 are 0.19% and 2.0%, respectively. An additional systematic error
results from the resources of the Analog Digital Conversion (ADC) that can be
approximated to be of the order of 0.15% for detector G resulting in a total
uncertainty of 0.25%. Thus variations on a 0.5% and 1% level have a
significance of 2 and 4$\sigma$, respectively. From the figure it is clear,
that single detectors with a count rate of $\approx$300,000 counts per 10
minutes can resolve cosmic ray variations that are as small as one percent in
less than an hour while the nominal coincidence channels (hundreds of counts
per minute) do not offer sufficient statistics for such an analysis. The
disadvantage of the single counter on the other side is the loss of any direct
information on the particle types and energies responsible for the count rate.
However, by calculating the energy dependent response for each single counter,
this information can be indirectly derived to some extend allowing the
counters for example to be used to study Galactic Cosmic Ray (GCR) modulation
(Kühl et al., 2015b). Figure 4 displays the daily averaged quiet time GCR flux
variation over the last two solar cycles and from July 2019 to July 2020,
respectively. The black and red line show the count rate variation of the anti
coincidence detector G and the SSD F, respectively. A count rate threshold in
SSD B was used to omit time periods of increased energetic particle fluxes
(e.g. SEPs). Although both curves follow each other during most times the
anti-coincidence detector G reflects some increases not measured by the SSD F.
The reason for this lies in the response of G to low energy electrons that are
not measured by SSD B. The red line reflects level of the highest galactic
cosmic ray flux ever reported during the space era (Zhao et al., 2014). From
the lower panel of the figure it is evident, that we reached a higher flux
level in 2020 leading to a harsher GCR caused radiation environment close to
Earth than in the previous solar minimum.

Figure 4: Uper and lower panel display the daily averaged quiet time galactic
cosmic ray flux variation over the last two solar cycles and from July 2019 to
July 2020, respectively. The black and red line show the count rate variation
of the anti coincidence detector G and the SSD F using counting rate
thresholds in SSD B omitting time periods of increased energetic particle
fluxes (for details see text).

In addition to particles, photons from solar flares can also cause energy
depositions in the detectors. Especially detector A with a threshold of only
30 keV is sensitive to X-ray. Figure 5 shows the count rate of the lower two
electron channels (red and teal) as well as the single counter rate of
detector A (blue) and B (green) during the SEP event on July 9, 1996 (see for
example Sierks et al., 1997). All channels and single counter show an increase
in count rate starting at 9:20 UT with a maximum around 9:40 UT followed by a
gradual decrease due to SEPs. However, the detector A single count rate (as
well as the lowest electron channel up to small extent) have another maximum
at approximately 9:10 UT lasting only for a couple of minutes that has been
caused by photons interacting in the detector (compare to Posner, 2007, and
references therein). The photons stem from an X2.6 x-ray flare located at S10
and W30 that had an onset, maximum and end time at 09:01 UT, 09:12 UT and
09:49 UT, respectively (from ftp://ftp.ngdc.noaa.gov/STP/space-weather/solar-
data/solar-features/solar-flares/x-rays/goes/xrs/). These capabilities of
detecting solar flare photons can be exploited for studies on and/or
forecasting schemes (see for example Núñez et al., 2018, and references
therein). It has to be noted, however, that the flare signal presented in
Figure 5 was measured during an early phase of the mission. Since the
background noise of detector A increased from roughly 1000 counts per minute
to the current level of 1-3 105 counts per minute over the first ten years of
the mission, the sensitivity of EPHIN to solar flare photons is restricted to
larger flares (M and X class) for the majority of the mission.

Figure 5: One minute resolution count rates of the electron channels E150
(red), E300 (teal), the single detectors A (blue) and B (green) on July 9,
1996.

####  2.1.2 Histograms

EPHIN provides histograms of the total energy deposition in the detector
stack. These histograms are given for four different particle penetration
depths (e.g. particles stopping in B, C, D or E) with 64 bins each
(corresponding to roughly 30 channels populated by proton and helium each in
total). Although the energy deposition is calculated using the PHA electronic
chain on-board, the histogram data is not to be confused with the PHA data
product which is an entirely independent data product. The histograms are
averaged over eight minutes due to bandwidth limitations. Furthermore, it is
not possible to distinguish between different isotopes (i.e. 3He/4He) based on
the histograms. However, the histogram data combines most advantages from PHA
data (higher energy resolution than coincidence channels, better particle
identification) with the good statistics from the counter data products (in
contrast to PHA, which only represents a statistical sample). Hence,
applications such as a velocity dispersion analysis (VDA) during the onset of
an SEP event (e.g. Krucker and Lin, 2000) can benefit from this data set.
While the time-resolution of 8 minutes is relatively large for a VDA, EPHINs
large energy range from 4.3 to 53 MeV allows to measure protons with
differences in arrival time of roughly 80 minutes (assuming a 1.3 AU parker
spiral) and hence, still provide meaningfull VDA results (for a VDA based on
EPHIN data see e.g. Klein, K.-L. and Posner, A., 2005). Furthermore, the
calculation of fluence spectra integrated over an entire SEP event (e.g.
Mewaldt et al., 2005) is possible based on this data set. While the Level 1
histograms only provide counts at the different energy deposition, a new data
product providing fluxes and sophisticated particle separation is under
preparation.

####  2.1.3 Coincidence Counter

The coincidence counters contain counting rates for individual particle
species plus an additional channel for penetrating particles (the integral
channel). For the four ranges stopping in the detector stack (particle
stopping either in B, C, D or E), different thresholds for detector A (given
in table 3, Müller-Mellin et al., 1995) are used to distinguish roughly
onboard between electrons, protons and helium particles (resulting in 12
stopping channels in total). Note that particles heavier than helium can not
be resolved by EPHIN due to the usage of linear amplifiers in the electronics
and hence, they will be counted in the helium channel. Since the contribution
of these heavier particles can be considered small with respect to the helium
flux (see Reames, 2019; George et al., 2009, for SEPs and GCRs, respectively)
the systematic uncertainty introduced by this false identification can be
neglected. The simple particle identification by the thresholds, however, has
significant flaws regarding the electron, proton and helium identification as
well. Electrons can spoil the proton channel, while the helium and electron
channels also have a finite response to protons. These small relative
contributions can get significant depending on the measured actual flux
environment (e.g. electron rich SEP events, or SEP events with a small helium
to proton ratio). The advantage of this data product is the higher statistics
in comparison to other data products with improved particle separation such as
the Level 3 data (section 2.3).
For the proton and helium measurements, the count rates are further divided
according to trajectory inclination relative to the sensor axis. This
information can be used for path-length corrections of the energy depositions
resulting in an improved particle separation (i.e. isotopic resolution).

####  2.1.4 Level 1 PHA data

The Level 1 Pulse Height Analysis (PHA) contains detailed information for a
statistical sample (a simple priority table is ensuring at least a couple of
helium particles) of the measured particles. For those particles covered, the
segments triggered in detectors A and B, the identified coincidence, as well
as the height of the voltage pulse measured by the EPHIN electronics that was
caused by the energy deposition of the particle and the used amplifier (low or
high gain) is given for each detector (for details see Müller-Mellin et al.,
1995). Hence, in contrast to the histogram data, not only the total energy
deposition of a particle in the detector stack but also the information
regarding the variations and changes in energy deposition along its trajectory
through the various detectors is available. However, we strongly recommend to
utilize the Level-2 PHA data sets as explained below.

###  2.2 Level 2

The Level 2 data consists of 1) flux data for electrons, protons and helium
particles with four energy channels each as well as one penetrating channel
(RL2) data as well as 2) PHA data in units of MeV. Both data sets can be
downloaded from http://ulysses.physik.uni-kiel.de/costep/level2/

####  2.2.1 RL2 data

Based on the Level 1 coincidence counters and simulations of the instrument
yielding the response function of the instrument for different particles and
energies, the flux is calculated for all twelve stopping channels (i.e. four
channels for electrons, protons and helium) in units of (cm2 sr s
MeV/nucleon)-1 plus the integral channel counting particles producing signals
in detectors A-F in units of (cm2 sr s)-1. While the application of such
response factors is straightforward for ions due to their almost constant
response over energy (see nominal proton response function in figure A1 in the
appendix from Kühl and Heber, 2019), the large amount of scattering of
electrons in the detector results in a more complex energy response typically
requiring a more detailed treatment such as the bow-tie analysis (Boudouridis
et al., 2020, and references therein). For the RL2 files, such an analysis was
not performed and the response has been only approximated by integration over
the energy dependent response functions resulting in systematic uncertainties
which are depending on the spectral shape of the observed electrons. Other
advantages and disadvantages of this data product include those given in
section 2.1.3.

####  2.2.2 Level 2 PHA data

Based on the Level 1 Pulse Height Analysis (PHA) and the calibration of the
instrument, the energy deposition in each detector can be calculated. The
Level 2 PHA contains this detailed information for a statistical sample (a
simple priority table is ensuring at least a couple of helium particles) of
the measured particles. For those particles covered, the segments triggered in
detectors A and B, the identified coincidence, as well as energy deposition in
units of MeV of the particle in each detectors given.
Utilizing the Level 2 PHA data, EPHIN is capable of providing energy spectra
with small & variable energy channels, identification and separation of
different isotopes (i.e. 3He/4He) as well as more complex, new data products
(see section 2.3 and 6).
As mentioned above the first two SSDs are sectorized in order to determine the
different path length of particles in order to improve the determination of
the 3He/4He ratio (Müller-Mellin et al., 1995). The PHA files store more
detailed information, i.e. which sector is hit in the first two SSDs for every
single particle evaluated. Banjac et al. (2015) has shown that 36 virtual
telescopes defined by the different sector combinations can be utilized to
investigate the particle pitch angle distribution in the field of view of the
instrument with an approach recently also proposed for the Parker Solar
Probe’s ISIS/EPI-Hi (McComas et al., 2016). Unfortunately, the counting
statistics of the PHA data only allows to determine the deviation from an
isotropic distribution for only a few events.

###  2.3 Level 3 Ion Fluxes

Due to increasing noise in detector D, EPHIN was commanded to Failure mode D
in October, 2017, (see section 1) limiting the energy resolution of the
nominal count rate based data products. Therefore, Kühl and Heber (2019) have
developed a method based on PHA data in order to restore the original energy
resolution for protons and helium using only the information of energy
depositions in detectors A, B and C. The data product has been successfully
validated against ACE/SIS, SOHO/ERNE and GOES11/EPS. However, due to the
methodology, neither electron fluxes nor isotopic resolution (i.e. 3He/4He)
could be retrieved by this method. Furthermore, the PHA based nature of the
data product can cause time-dependent variations of the statistical
uncertainty (i.e. at the onset of an SEP event, when electrons are filling
most of the PHA buffer). Nonetheless, the long-term availability of the
measurements make the data especially interesting for space weather studies
and monitoring (cf. Kühl and Heber, 2019). The resulting data product offers
consistent measurements from 1995 on, various time resolutions (from 1 minute
to daily averages) as well as statistical and systematical uncertainty in a
user friendly file format accessible here: http://ulysses.physik.uni-
kiel.de/costep/level3/l3i/. Furthermore, the method has been extended such
that the flux of protons and helium nuclei can be derived for any energy
channel in the range from 4 to 53 MeV/nucleon allowing an easy inter-
calibration between various different missions.

###  2.4 Penetrating particles

Figure 6: Single Event upsets recorded by SOHO’s solid state recorder and
EPHIN measured flux of 500-700 MeV protons.

Kühl et al. (2015a) have presented a method in order to extend the energy
range of EPHIN for protons up to energies of above 1 GeV based on the Level 2
PHA data. Utilizing extensive GEometry ANd Tracking (GEANT4)
(Geant4-Collaboration, 2003) simulations of the instrument, it has been shown
that the initial energy of a given proton penetrating the instrument (i.e.
triggering detectors A to F) can be derived from the energy losses in SSD C
and D. Due to too high noise in detector D (ultimately leading to the switch
on of the EPHIN failure mode D), the method is restricted to data before 2017.
However, those high energy measurements of GCRs (Kühl et al., 2016) as well as
SEPs (Kühl et al., 2017) have been proven to be not only of scientific value
for the study of GCR modulation as well as SEP acceleration and transport, but
can be also utilized for space weather applications. As an example, Figure 6
presents the flux of 500-700 MeV protons measured by EPHIN in comparison to
the single event upsets (SEUs) recorded by SOHO’s solid state recorder (see
Curdt and Fleck, 2015, for details). The dependency of the SEU rate on the
solar cycle variations due to the GCR modulation as well as sporadic increases
due to SEP events can be clearly seen in this comparison. Annual averaged
proton flux spectra from 250 MeV to 1.6 GeV based on EPHIN measurements can be
downloaded from the cosmic ray database (Maurin et al., 2014):
https://lpsc.in2p3.fr/cosmic-rays-db/.

##  3 Space Weather Applications

###  3.1 Long-term Data / Space Weather Monitoring

EPHIN data has been used in various SEP event analysis (e.g. Cohen and
Mewaldt, 2018; Dresing et al., 2018; Gómez-Herrero et al., 2002; Lario et al.,
2014) and has contributed data to multiple SEP event catalogues (Vainio et
al., 2013; Paassilta et al., 2017).
Furthermore, the availability of EPHIN data over more than two solar cycles
also allows for long-term studies of the variation of the GCR flux over two
solar cycles of different magnetic polarity with special emphasis on the
contribution of SEP events to the overall fluence observed during this long-
term mission. Kühl and Heber (2019) has shown that for protons and helium
particles the fluence during as little as 10 days contributes between 24% and
79% to the overall fluence during 23 years of measurements depending on the
particles energy. The large SEP events responsible for these contributions
have been also linked to significant solar array degradation and SEUs in
SOHO’s electronics (Curdt and Fleck, 2015).
The derived flux for high energetic penetrating particles has been also used
for validation purposes of the GLE forecast scheme HESPERIA UMASEP-500 (Núñez
et al., 2017).

###  3.2 Near Real Time Data / Space Weather Forecast

While most EPHIN data sets are only available with a couple of days of delay,
some near real time Level 1 data are available and includes SCI and PHA files
for each available minute individually (in contrast to the daily merged files
of all other data products). However, due to limited downlink possibilities,
these near real time data are only available for a couple of hours each day
depending on the ground station coverage (Malandraki and Crosby, 2018) . The
data is provided here: http://ulysses.physik.uni-
kiel.de/costep_realtime/realtime/lvl1/

Figure 7: Performance of the HESPERIA REleASE forecast scheme based on non-
realtime data from EPHIN (left) and ACE/EPAM (right) for GLE 71 in May, 2012.

One of the space weather related applications of the near real time data of
EPHIN is the forecasting scheme Relativistic Electron Alert System for
Exploration (REleASE) introduced by Posner (2007). SEP events can cause an
increase of the particle flux by several orders of magnitude for various
particle types such as electron, protons and ions. While the radiation hazard
caused by these particle populations has a strong dependency on the event
characteristics such as the energy spectra and abundances of the various
particle types as well as the location of the observer and the surrounding
shielding, it can be generally assumed that the relativistic, minimum-ionizing
electrons have a significant lower impact than the protons and ions. However,
since electrons and protons are often accelerated and injected into the
interplanetary medium at the same time during an SEP event, the shorter
propagation time of the electrons can be exploited in order to forecast the
proton that arrive at L1 with an delay of approximately 30-90 minutes
depending on the interplanetary propagation as well as the particle energy. In
detail, the REleASE system utilized an empirical relationship between the
expected proton flux at energies $>$30-50 MeV measured 30, 60 or 90 minutes
later than the electrons dependent on measured electron flux and its increase
that have been derived based on historical EPHIN data. The derived forecasting
matrix can then be applied to near real time electron measurements in order to
issue warnings and forecasts of the proton flux at 1AU. Posner and Strauss
(2020) have also studied the feasibility of a REleASE system supplemented with
data from a second instrument near Mars in order to provide warnings for a
possible future human mission to Mars.
The REleASE scheme has been recently further improved as part of the EU
HORIZON 2020 HESPERIA (High Energy Solar Particle Events forecasting and
Analysis) project (Malandraki and Crosby, 2018). Besides software updates and
performance studies based on historical data that have not been used in the
calculation of the forecasting matrix, the HESPERIA REleASE scheme has been
developed such that its innovation is to enable derivation and provision of
the 30-50 MeV proton forecast based on electron measurements from ACE/EPAM
(Gold et al., 1998) as well as EPHIN. The advantage of ACE/EPAM is the higher
time coverage of the near real time data. HESPERIA REleASE provides forecasts
in a complementary way based on both instruments in parallel. As an example of
the forecast performance, Figure 7 shows the May, 2012 SEP event and the
related forecast based on historical data for EPHIN (left) and EPAM (right).
The upper panel presents the proton flux measured by EPHIN while the three
central panels show the proton forecast for 30, 60 and 90 minutes and the
bottom panels shows the electron flux used for the forecast. From the figure
it is clear that HESPERIA REleASE based on either of the two instruments was
able to forecast the event.
The real time plot as well as a collection of historical plots can be found on
the HESPERIA website:
https://www.hesperia.astro.noa.gr/index.php/results/real-time-prediction-
tools/release
It is noteworthy that the HESPERIA REleASE model has been selected on
international level and participates in the Integrated Solar Energetic Proton
Alert/Warning System (ISEP) (https://ccmc.gsfc.nasa.gov/isep/). ISEP is a
3-year collaboration project between the Space Radiation Analysis Group (SRAG)
at Johnson Space Center and the Community Coordinated Modeling Center (CCMC)
at Goddard Space Flight Center to bring state-of-the-art space weather models
from research and development at universities and small businesses to
operational use at NASA (R2O). These models have a user interface in the form
of the SEP scoreboard (https://ccmc.gsfc.nasa.gov/challenges/sep.php) that
allows the SRAG console operator to view and compare the results from several
different models simultaneously. HESPERIA REleASE has been recently
successfully integrated in the CCMC/SEP scoreboard, providing unique space
weather predictions for crew protection in human spaceflight in the exo-LEO
mission era.

##  4 Lessons Learned & Possible Improvements

Figure 8: Solar energetic particle event onset, observed over 12 hours on Nov.
4, 2001 by GOES and SOHO. The top panel shows the GOES/SXT low-energy X-ray
flux. The second panel from top shows a spectrogram of energetic electrons in
the 0.25-10 MeV range observed by SOHO/COSTEP-EPHIN. Below that is displayed a
4-50 MeV proton spectrograms observed by EPHIN, and at the bottom the SEPEM
energy channels of GOES-8/SEM displayed as a spectrogram. Only energy channels
with full overlap between GOES 8/SEM and EPHIN are shown.

The well thought-out and robust instrument design as well as the carefully
crafted data products have resulted in EPHIN being a valuable asset for
particle detection not only for the space science but also for the space
weather community. Its large geometric factor and low background level
compared to other instruments (e.g. GOES/EPS, see fig. 5 from Kühl and Heber,
2019) allow for studies of even small SEP events, additional data products
such as the single counters have further extended the measurement capabilities
and more than a billion particles evaluated and stored as PHA data allow for
detailed studies over the 25 years of measurements. However, there are also
some lessons learned and possible improvements for EPHIN that should be taken
into account during the design of the next generation of particle instruments.
The extreme mixed radiation field of relativistic electrons, protons and heavy
ions, X- and gamma rays, and secondary particles generated from the
interaction of primary particles with matter pose significant challenges to
energetic particle detectors such as EPHIN. This problem is illustrated in
Figure 8, which shows spectrograms of energetic particles measured with
SOHO/COSTEP-EPHIN and GOES-8/SEM, and GOES X ray observations originating from
the associated active region flare. Issues affecting EPHIN and SEM have
already described for this event in the Appendix of Posner (2007). SEM
essentially measures integral fluxes of energetic particles in a broad energy
range, spanning from a few MeV to almost 1 GeV. It uses passive shielding and
a crude way of identifying counts for the various integral channel (such as
$>$10 MeV). Unfortunately, the passive shielding does not stop protons
exceeding 60 MeV from entering the detector. Multiple teams have attempted
cleaning up the GOES measurements over time. Zwickl (unpublished; see Appendix
A. in Rodriguez et al., 2017, for details) created differential channels by
using a subtraction scheme of different integral channels, but admitted that
difficulties arise that the 15–40 MeV proton channel sees with a large
geometric factor ions that penetrate its passive shielding. The Zwickl-derived
GOES 8/SEM for this event are shown in Figure 1 of Posner (2007). As compared
to the cleaner EPHIN data set, which clearly shows that the proton event
starts with an energy dispersion signature, the GOES-8/SEM onset in the 15-40
MeV channel is off by 30 min.
A more recent development is the SEPEM (sepem.eu, Jiggens et al., 2012;
Sandberg et al., 2014) treatment of GOES energetic proton fluxes. The output
of SEPEM for the Nov. 4, 2001 event is shown in the bottom panel of Figure 8
as a spectrogram consisting of all differential proton channels that have full
overlap with the EPHIN proton energy range. It is apparent that the SEPEM-
generated onset of this event contains a large amount of high-energy protons
that spill over into the 5-40 MeV proton channels, completely obscuring the
actual onset of the new event. Moreover, the SEPEM spectra during this time
frame suggest a decreasing intensity with energy in this energy range, which
would contradict the cleaner measurement made by EPHIN. The superior
performance of EPHIN during this critical onset time is thanks to the active
anti-coincidence system, an essential component of any accurate SEP detection
system that is supposed to provide reliable measurements while being exposed
to any given particle environment. At the same time, the suppression of pre-
event counts is due to the mixed radiation field EPHIN is exposed to,
particularly the onset of the intense electron event. As mentioned in Section
2.1.3, the threshold of the coincidence system responsible for the initial
particle identification scheme is set too low, allowing a large number of
electrons being recorded in the low-energy proton channels. This combined with
the low count rate of actual protons from the previous SEP event suppresses
the number of pulse-height analyzed protons in this channel and replaces them
with electrons, which are sorted out by the subsequent dE/dx vs E analysis.
Only very few actual proton counts remain until the new SEP proton event kicks
in, generating the well-known energy dispersion pattern of fast protons
arriving before lower-energy slower protons. Any vertical column in the EPHIN
spectrogram is representative of a momentary (2-min) spectrum of the event. It
is clearly seen that the spectra are anomalous, i.e., the spectral intensity
increases with energy below the dispersion onset, and above that rolls over
into a nominal proton spectrum that decreases in intensity with energy. In a
new detector system, this effect can effectively be reduced by adjusting the
lower-energy proton channel threshold or by utilizing increased on-board
processing capabilities that allow for more complicated particle separation
rather than just one threshold of one detector, which would significantly
reduce the spill-over of electrons into the proton channels. None of the
energy dispersion phenomena are captured in the GOES 8/SEM data set shown here
as part of the SEPEM output.
While the dynamic range of EPHIN achieved with the automatic switch-off of the
outer segments of detectors A and B during periods of high fluxes (resulting
in a decrease of the geometric factor by a factor $\approx$ 24) is superior
compared to other instruments (e.g. SOHO/ERNE, see fig. 7 from Kühl and Heber,
2019), it does not prevent EPHIN from suffering pile-up and dead-time issues
during the largest SEP events. Hence, a reduction of the size of the inner
segments and implementing a sophisticated livetime counter in the electronics
of EPHIN would have improved the data quality during these events.
Furthermore, additional material surrounding the anti-coincidence detector G
would have decreased the high count rates of this detector dominated by low
energy particles during SEP events and thus further reducing the dead-time of
the instrument. It has to be also noted, that the ring switch-off can
introduce additional uncertainties regarding the particle identification.
While the ring is switched off, particles penetrating an outer segment of
detector A and hitting the central segment of B at the same time as a valid
particle hitting both central segments will add its energy losses in B and
deeper detectors to the measurement, artificially increasing the measured
energy of the detected particle and therefore resulting in a too hard
spectrum. This effect could explain the hard electron spectrum in figure 8,
while the sudden decrease in flux above 600 keV is most likely due to an
underestimation of the dead layer of detector C. Therefore, it should also be
considered for future instruments to add a complimentary high flux instrument,
also based on EPHIN principles, with a lower geometric factor and sufficient
passive and active anti-coincidence to allow measurements of particle fluxes
during the most extreme SEPs events. While this second instrument would not be
able to detect small SEP events, the combination of them would increase the
dynamic range of the observable fluxes. The realization of such a design can
be inspired by the pixel detector used for the PSP/ISIS/EPI-Hi (McComas et
al., 2016). The possibility of using two or more instruments with adjustable
geometric factor similar to EPHIN, each optimized for measuring only in a
small energy range should be also explored as an alternative approach.
GEANT4 simulations show that $\delta$-electrons generated by ions stopping in
detector A can propagate to detector B and cause energy depositions above the
threshold of 60 keV. This issue of particles below the nominal energy range
being counted by the coincidence logic could have been easily resolved by
increasing the threshold to above 100 keV.
It is moreover recognized that proton energies of 30-100 MeV are of
significant relevance for radiation safety, as they penetrate deeper in matter
and thus can more readily reach the inside of exploration vehicles as well as
astronauts in extravehicular suits than protons of lower energies. The EPHIN
limitation of stopping protons of $<$54 MeV could have been mitigated by
either using more and/or thicker detectors at the rear of the particle
telescope or by including a scintillation detector in the detector stack. The
former solution provides better energy resolution on the costs of requiring
more channels or high bias voltage for thicker detectors. These thicker
detectors, such as the lithium drifted ones used for EPHIN, are also less
reliable for a long-term mission (see section 1). McComas et al. (2016),
however, have shown that two successive detectors connected to the same front-
end electronics as used for the PSP/ISIS/EPI-Hi instrument can mitigate these
downsides. Utilizing a scintillation detector as used for i.e. the Energetic
Particle Detector on Solar Orbiter (Rodriguez-Pacheco, J. et al., 2019)
provides a high stopping power especially for heavier ions without the need of
several channels while reducing the energy resolution and introducing
systematic uncertainties due to higher temperature dependencies as well as
non-linearities such as quenching.
Histograms that store the count rate distribution for protons, Helium but also
for energetic electrons in all 36 virtual telescopes (cf. section 2.2.2) would
have allowed to compute the directional distribution without the statistical
limitation of an PHA based analysis.
The calculation of the high energy spectra of penetrating particles is limited
by statistics since it is based on PHA data. This limitation could have been
eliminated by providing an energy loss histogram for particles penetrating the
entire detector stack in the histogram data. However, the method inhabits a
significant systematic uncertainty due to EPHIN’s inability of separating
forward and backward penetrating particles as well as relativistic protons and
electrons (Kühl et al., 2016). The E6 Instruments on the Helios spacecrafts
have shown that adding a Cherenkov detector behind the detector stack can be
utilized in order to distinguish between both directions (Marquardt and Heber,
2019) while the COSPIN/KET instrument on Ulysses (Simpson et al., 1992) has
proven the possibility to distinguish between relativistic protons up to 2.1
GeV and electrons utilizing an Aerogel Cherenkov detector.
An energy-loss histogram of one of the inner detectors in anti-coincidence to
all other detectors would have provided more capabilities in the analysis of
$\gamma$-rays following the approach of the Lunar Lander Neutron and Dosimetry
experiment (LND) on Chang’e4 (Wimmer-Schweingruber et al., 2020) and
PSP/ISIS/EPI-Hi (McComas et al., 2016).
By utilizing the increased on-board processing capabilities of these state of
the art instruments compared to the 25 year old EPHIN electronics, the
proposed improvements could be included without the need of significant higher
telemetry rate. For example, the 2.8 Gbit per orbit telemetry from
PSP/ISIS/EPI-Hi translates to an average of 192 bit per second for the first
orbit of 168 days which is similar to EPHINs nominal telemetry rate of 172 bit
per second. The mass and power budget of such an improved instrument, however,
would depend on the detailed design and is therefore difficult to estimate.

##  5 Summary and Outlook

Providing 25 years of continuous measurements during two entire solar cycles
and a variety of different solar activities resulting in different SEP
environments, EPHIN on-board SOHO is a valuable asset regarding particle
detection for the space weather and solar science communities. Due to the low
background thanks to the lateral anti-coincidence and the large geometric
factor as well as the availability of PHA data (including penetrating
particles), it was possible to provide new data products with enhanced
scientific value utilizing revised simulations (GEANT4) of the already 25 year
old instrument. The numerous data products summarized here have been used for
multiple long-term cosmic ray and SEP studies as well as space weather
applications such as SEP event forecast and the analysis of single event
upsets in spacecraft electronics.
On the technical level, a number of benefits and limitations of EPHIN have
been pointed out, which a future instrument best qualified for space weather
purposes should incorporate and overcome, respectively.
In order to accurately evaluate SEP forecasting systems, it is imperative that
the actual onsets of SEP events can be determined, both in terms of actual
event onset by differential proton channel, and in terms of when a
differential proton channel exceeds a threshold level that is relevant for
protecting human explorers in deep space. EPHIN is currently one of only few
detectors that have this capability due to the large dynamic range (which
could be further improved by adding a second high flux instrument) as well as
the low background and spill-over of different particle types mainly
accomplished by the anticoincidence system. Furthermore, a variety of data
products that have extended the measurement capabilities of EPHIN have been
derived and provided to the community based on large amount of PHA data.
Improvements would include a slightly higher energy range (up to $\approx$100
MeV) achieved by adding more SSDs, a complex and sophisticated onboard
particle separation as well as further passive material surrounding the active
parts in order to decrease pile up effects and deadtime limitations. The
possibility of an enhanced detection of neutral particles should also be taken
into account.
We propose a Proton Electron Comprehensive Telescope for Moon/Mars Exploration
(ProtECT-ME) based on the EPHIN design that would incorporate these features
and improvements in order to be best suited for radiation monitoring and
providing real time data for SEP forecast schemes.
With further improvements such as a rejection system for backward penetrating
particles and an improved separation between relativistic electrons and
protons as well as onboard data products such as histograms for penetrating
particles, such an instrument could be also used for aviation SEP warnings.

###### Acknowledgements.

The SOHO/EPHIN project is supported under grant 50 OC 1702 by the German
Bundesministerium für Wirtschaft through the Deutsches Zentrum für Luft- und
Raumfahrt (DLR). This project has received funding from the European Union’s
Horizon 2020 research and innovation program under grant agreement No. 637324.
RGH acknowledges the financial support by the Spanish Ministerio de Ciencia,
Innovación y Universidades FEDER/MCIU/AEI Projects ESP2017-88436-R and
PID2019-104863RB-I00.

## References

  * Banjac et al. (2015) Banjac, S., R. Gómez-Herrero, B. Heber, P. Kühl, and C. Terasa, 2015.  EPHIN anisotropy measurement capability.  In Journal of Physics Conference Series, vol. 632 of _Journal of Physics Conference Series_ , 012048.  10.1088/1742-6596/632/1/012048.
  * Boudouridis et al. (2020) Boudouridis, A., J. V. Rodriguez, B. T. Kress, B. K. Dichter, and T. G. Onsager, 2020.  Development of a Bowtie Inversion Technique for Real-Time Processing of the GOES-16/-17 SEISS MPS-HI Electron Channels.  _Space Weather_ , 18(4), e02403.  10.1029/2019SW002403.
  * Cohen and Mewaldt (2018) Cohen, C. M. S., and R. A. Mewaldt, 2018.  The Ground-Level Enhancement Event of September 2017 and Other Large Solar Energetic Particle Events of Cycle 24.  _Space Weather_ , 16(10), 1616–1623.  10.1029/2018SW002006.
  * Curdt and Fleck (2015) Curdt, W., and B. Fleck, 2015.  Solar and Galactic Cosmic Rays Observed by SOHO.  _Central European Astrophysical Bulletin_ , 39, 109–118.  1505.07952.
  * Domingo et al. (1995) Domingo, V., B. Fleck, and A. I. Poland, 1995.  The SOHO Mission: an Overview.  _Sol. Phys._ , 162(1-2), 1–37.  10.1007/BF00733425.
  * Dresing et al. (2018) Dresing, N., R. Gómez-Herrero, B. Heber, A. Klassen, M. Temmer, and A. Veronig, 2018.  Long-lasting injection of solar energetic electrons into the heliosphere.  _A &A_, 613, A21.  10.1051/0004-6361/201731573, 1802.04722.
  * Geant4-Collaboration (2003) Geant4-Collaboration, 2003.  Geant4 - a simulation toolkit.  _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , 506(3), 250 – 303.  10.1016/S0168-9002(03)01368-8, URL http://dx.doi.org/10.1016/S0168-9002(03)01368-8.
  * George et al. (2009) George, J. S., K. A. Lave, M. E. Wiedenbeck, W. R. Binns, A. C. Cummings, et al., 2009.  Elemental Composition and Energy Spectra of Galactic Cosmic Rays During Solar Cycle 23.  _ApJ_ , 698(2), 1666–1681.  10.1088/0004-637X/698/2/1666.
  * Gold et al. (1998) Gold, R. E., S. M. Krimigis, I. Hawkins, S. E., D. K. Haggerty, D. A. Lohr, E. Fiore, T. P. Armstrong, G. Holland, and L. J. Lanzerotti, 1998\.  Electron, Proton, and Alpha Monitor on the Advanced Composition Explorer spacecraft.  _Space Sci. Rev._ , 86, 541–562.  10.1023/A:1005088115759.
  * Gómez-Herrero et al. (2002) Gómez-Herrero, R., M. D. Rodríguez-Frías, L. del Peral, R. Müller-Mellin, and H. Kunow, 2002.  SOHO/EPHIN observation of a multiple large solar energetic particles event in November 1997.  _Astroparticle Physics_ , 17(1), 1–12.  10.1016/S0927-6505(01)00136-0.
  * Heber et al. (2015) Heber, B., C. Wallmann, D. Galsdorf, Herbst K., P. Kühl, et al., 2015\.  Forbush decreases associated to Stealth Coronal Mass Ejections.  _Central European Astrophysical Bulletin_ , 39, 75–82.
  * Howard and Harrison (2013) Howard, T. A., and R. A. Harrison, 2013.  Stealth Coronal Mass Ejections: A Perspective.  _Sol. Phys._ , 285(1-2), 269–280.  10.1007/s11207-012-0217-0.
  * Jiggens et al. (2012) Jiggens, P. T. A., S. B. Gabriel, D. Heynderickx, N. Crosby, A. Glover, and A. Hilgers, 2012.  ESA SEPEM Project: Peak Flux and Fluence Model.  _IEEE Transactions on Nuclear Science_ , 59(4), 1066–1077.  10.1109/TNS.2012.2198242.
  * Klein, K.-L. and Posner, A. (2005) Klein, K.-L., and Posner, A., 2005.  The onset of solar energetic particle events: prompt release of deka-MeV protons and associated coronal activity.  _A &A_, 438(3), 1029–1042.  10.1051/0004-6361:20042607, URL https://doi.org/10.1051/0004-6361:20042607.
  * Krucker and Lin (2000) Krucker, S., and R. P. Lin, 2000.  Two Classes of Solar Proton Events Derived from Onset Time Analysis.  _ApJ_ , 542(1), L61–L64.  10.1086/312922.
  * Kühl et al. (2015a) Kühl, P., S. Banjac, N. Dresing, R. Goméz-Herrero, B. Heber, A. Klassen, and C. Terasa, 2015a.  Proton intensity spectra during the solar energetic particle events of May 17, 2012 and January 6, 2014.  _A &A_, 576, A120.  10.1051/0004-6361/201424874.
  * Kühl et al. (2015b) Kühl, P., S. Banjac, B. Heber, J. Labrenz, R. Müller-Mellin, and C. Terasa, 2015b.  Extended Measurement Capabilities of the Electron Proton Helium INstrument aboard SOHO - Understanding single detector count rates.  _Central European Astrophysical Bulletin_ , 39, 119–124.
  * Kühl et al. (2017) Kühl, P., N. Dresing, B. Heber, and A. Klassen, 2017.  Solar Energetic Particle Events with Protons Above 500 MeV Between 1995 and 2015 Measured with SOHO/EPHIN.  _Sol. Phys._ , 292(1), 10.  10.1007/s11207-016-1033-8, 1611.03289.
  * Kühl et al. (2016) Kühl, P., R. Gómez-Herrero, and B. Heber, 2016.  Annual Cosmic Ray Spectra from 250 MeV up to 1.6 GeV from 1995 - 2014 Measured with the Electron Proton Helium Instrument onboard SOHO.  _Sol. Phys._ , 291(3), 965–974.  10.1007/s11207-016-0879-0, 1603.00676.
  * Kühl and Heber (2019) Kühl, P., and B. Heber, 2019.  Revising More Than 20 Years of EPHIN Ion Flux Data—A New Data Product for Space Weather Applications.  _Space Weather_ , 17(1), 84–98.  10.1029/2018SW002114.
  * Lario et al. (2014) Lario, D., N. E. Raouafi, R. Y. Kwon, J. Zhang, R. Gómez-Herrero, N. Dresing, and P. Riley, 2014.  The Solar Energetic Particle Event on 2013 April 11: An Investigation of its Solar Origin and Longitudinal Spread.  _ApJ_ , 797(1), 8.  10.1088/0004-637X/797/1/8, 1410.5490.
  * Malandraki and Crosby (2018) Malandraki, O. E., and N. B. Crosby, 2018.  Solar Particle Radiation Storms Forecasting and Analysis, vol. 444.  Springer.  10.1007/978-3-319-60051-2.
  * Marquardt and Heber (2019) Marquardt, J., and B. Heber, 2019.  Galactic cosmic ray hydrogen spectra and radial gradients in the inner heliosphere measured by the HELIOS Experiment 6.  _A &A_, 625, A153.  10.1051/0004-6361/201935413, 1905.01052.
  * Maurin et al. (2014) Maurin, D., F. Melot, and R. Taillet, 2014.  A database of charged cosmic rays.  _A &A_, 569, A32.  10.1051/0004-6361/201321344, 1302.5525.
  * McComas et al. (2016) McComas, D. J., N. Alexander, N. Angold, S. Bale, C. Beebe, et al., 2016\.  Integrated Science Investigation of the Sun (ISIS): Design of the Energetic Particle Investigation.  _Space Sci. Rev._ , 204(1-4), 187–256.  10.1007/s11214-014-0059-1.
  * Mewaldt et al. (2005) Mewaldt, R. A., C. M. S. Cohen, A. W. Labrador, R. A. Leske, G. M. Mason, M. I. Desai, M. D. Looper, J. E. Mazur, R. S. Selesnick, and D. K. Haggerty, 2005.  Proton, helium, and electron spectra during the large solar particle events of October-November 2003.  _Journal of Geophysical Research (Space Physics)_ , 110(A9), A09S18.  10.1029/2005JA011038.
  * Müller-Mellin et al. (1995) Müller-Mellin, R., H. Kunow, V. Fleißner, E. Pehlke, E. Rode, et al., 1995\.  COSTEP - Comprehensive Suprathermal and Energetic Particle Analyser.  _Solar Physics_ , 162(1), 483–504.
  * Núñez et al. (2018) Núñez, M., K.-L. Klein, B. Heber, O. E. Malandraki, P. Zucca, J. Labrens, P. Reyes-Santiago, P. Kuehl, and E. Pavlos.  HESPERIA Forecasting Tools: Real-Time and Post-Event, vol. 444 of _Astrophysics and Space Science Library_ , 113–131.  Springer, 2018.  10.1007/978-3-319-60051-2_7.
  * Núñez et al. (2017) Núñez, M., P. J. Reyes-Santiago, and O. E. Malandraki, 2017.  Real-time prediction of the occurrence of GLE events.  _Space Weather_ , 15(7), 861–873.  10.1002/2017SW001605.
  * Paassilta et al. (2017) Paassilta, M., O. Raukunen, R. Vainio, E. Valtonen, A. Papaioannou, et al., 2017.  Catalogue of 55-80 MeV solar proton events extending through solar cycles 23 and 24.  _Journal of Space Weather and Space Climate_ , 7, A14.  10.1051/swsc/2017013, 1707.00498.
  * Policella et al. (2009) Policella, N., H. Oliveira, and T. Siili, 2009.  SKEYP: AI Applied to SOHO Keyhole Operations.  In 2009 Third IEEE International Conference on Space Mission Challenges for Information Technology, 58–65.
  * Posner (2007) Posner, A., 2007.  Up to 1-hour forecasting of radiation hazards from solar energetic ion events with relativistic electrons.  _Space Weather_ , 5(5), 05001.  10.1029/2006SW000268.
  * Posner and Strauss (2020) Posner, A., and R. Strauss, 2020.  Warning Time Analysis From SEP Simulations of a Two-Tier REleASE System Applied to Mars Exploration.  _Space Weather_ , 18(4), e2019SW002,354.  E2019SW002354 2019SW002354, 10.1029/2019SW002354, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019SW002354, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019SW002354.
  * Reames (2019) Reames, 2019.  Element Abundances of Solar Energetic Particles and the Photosphere, the Corona, and the Solar Wind.  _Atoms_ , 7(4), 104.  10.3390/atoms7040104, 1910.01209.
  * Richardson et al. (1996) Richardson, I. G., G. Wibberenz, and H. V. Cane, 1996.  The relationship between recurring cosmic ray depressions and corotating solar wind streams at <=1 AU: IMP 8 and Helios 1 and 2 anticoincidence guard rate observations.  _J. Geophys. Res._ , 101(A6), 13,483–13,496.  10.1029/96JA00547.
  * Rodriguez et al. (2017) Rodriguez, J. V., I. Sandberg, R. A. Mewaldt, I. A. Daglis, and P. Jiggens, 2017.  Validation of the effect of cross-calibrated GOES solar proton effective energies on derived integral fluxes by comparison with STEREO observations.  _Space Weather_ , 15(2), 290–309.  10.1002/2016SW001533.
  * Rodriguez-Pacheco, J. et al. (2019) Rodriguez-Pacheco, J., Wimmer-Schweingruber, R. F., Mason, G. M., Ho, G. C., Sánchez-Prieto, S., Prieto, M., and et al., 2019.  The Energetic Particle Detector. Energetic particle instrument suite for the Solar Orbiter mission.  _A &A_.  10.1051/0004-6361/201935287, URL https://doi.org/10.1051/0004-6361/201935287.
  * Sandberg et al. (2014) Sandberg, I., P. Jiggens, D. Heynderickx, and I. A. Daglis, 2014.  Cross calibration of NOAA GOES solar proton detectors using corrected NASA IMP-8/GME data.  _Geophys. Res. Lett._ , 41(13), 4435–4441.  10.1002/2014GL060469.
  * Sierks et al. (1997) Sierks, H., I. Elent, W. Dröge, G. Mann, H. Aurass, H. Kunow, R. Müller-Mellin, J. Sequeiros, and S. McKenna-Lawlor, 1997.  Electron Spectra and Associated Type II and Type III Radio Bursts of the July 9, 1996 Solar Event.  In International Cosmic Ray Conference, vol. 1 of _International Cosmic Ray Conference_ , 297.
  * Simpson et al. (1992) Simpson, J. A., J. D. Anglin, A. Balogh, M. Bercovitch, J. M. Bouman, et al., 1992.  The ULYSSES Cosmic Ray and Solar Particle Investigation.  _A &AS_, 92(2), 365–399.
  * Vainio et al. (2013) Vainio, R., E. Valtonen, B. Heber, O. E. Maland raki, A. Papaioannou, et al., 2013.  The first SEPServer event catalogue ~68-MeV solar proton events observed at 1 AU in 1996-2010.  _Journal of Space Weather and Space Climate_ , 3, A12.  10.1051/swsc/2013030.
  * Vandenbussche (2006) Vandenbussche, F. C., 2006.  Saving SOHO.  In SOHO-17. 10 Years of SOHO and Beyond, vol. 617 of _ESA Special Publication_ , 29.
  * Wimmer-Schweingruber et al. (2020) Wimmer-Schweingruber, R. F., J. Yu, S. I. Böttcher, S. Zhang, S. Burmeister, et al., 2020.  The Lunar Lander Neutron and Dosimetry (LND) Experiment on Chang’E 4\.  _arXiv e-prints_ , arXiv:2001.11028.  2001.11028.
  * Zhao et al. (2014) Zhao, L. L., G. Qin, M. Zhang, and B. Heber, 2014.  Modulation of galactic cosmic rays during the unusual solar minimum between cycles 23 and 24.  _Journal of Geophysical Research (Space Physics)_ , 119(3), 1493–1506.  10.1002/2013JA019550, 1310.7076.

  *[CME]: Coronal Mass Ejection
  *[SEP]: Solar energetic particle
  *[SOHO]: SOlar and Heliospheric Observatory
  *[EPHIN]: Electron Proton Helium Instrument
  *[SSDs]: Solid State Detector
  *[PHA]: Pulse Height Analysis
  *[SSD]: Solid State Detector
  *[CMEs]: Coronal Mass Ejection
  *[FDs]: Forbush Decrease
  *[FD]: Forbush Decrease
  *[ADC]: Analog Digital Conversion
  *[GCR]: Galactic Cosmic Ray
  *[SEPs]: Solar energetic particle
  *[VDA]: velocity dispersion analysis
  *[GCRs]: Galactic Cosmic Ray
  *[GEANT4]: GEometry ANd Tracking


# A sparse approximate inverse for triangular matrices based on Jacobi
iteration

Zhongjie Lu Dipartimento di Matematica, Università di Pavia, via Ferrata 5,
27100, Pavia, Italy. Current address: School of Mathematical Sciences,
University of Science and Technology of China, Hefei, Anhui, 230026, China.
Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a simple sparse approximate inverse for triangular
matrices (SAIT). Using the Jacobi iteration method, we obtain an expression of
the exact inverse of triangular matrix, which is a finite series. The SAIT is
constructed based on this series. We apply the SAIT matrices to iterative
methods with ILU preconditioners. The two triangular solvers in the ILU
preconditioning procedure are replaced by two matrix-vector multiplications,
which can be fine-grained parallelized. We test this method by solving some
linear systems and eigenvalue problems with preconditioned iterative methods.

Keywords: triangular solver, ILU factorization, iterative method, linear
system, eigenvalue problem

##  1 Introduction

The incomplete LU (ILU) factorization is a type of general-purpose
preconditioning techniques for sparse linear systems. There are two main
problems in its parallelization. The first is the parallel generation of ILU
factors. Many theories and techniques have been used to improve its
parallelization. Here are some references on this problem [1, 3, 9, 14, 15,
18]. The second problem is solving triangular systems in the preconditioning
procedure. It is easy to solve exactly a triangular system using forward or
backward substitution method. However, this is a highly sequential process,
and it may be executed many times in solving a system. This is an obstacle to
exploit the performance of a parallel computing platform sufficiently. This is
the main problem we shall study in this paper.

The simplest parallelization strategy in solving triangular systems is to
compute the summations in the substitution methods in parallel:

$\begin{split}x_{i}=\frac{1}{a_{ii}}\left(b_{i}-\sum_{k=1}^{i-1}a_{ik}x_{k}\right)\quad\text{or}\quad
x_{i}=\frac{1}{a_{ii}}\left(b_{i}-\sum_{k=i+1}^{n}a_{ik}x_{k}\right).\end{split}$

This is a low-level concurrent method, especially for sparse triangular
matrices. Its corresponding blocking method can improve the performance. The
level-scheduling method with reordering techniques improve the parallelism
further [2]. These methods are exact solvers for triangular systems.

Another type of strategy is to replace the exact solutions of triangular
systems by approximate ones. This ideal is based on the fact that the
incomplete LU factorization is incomplete. This means that even the exact
solutions of the triangular systems in ILU are not ’exact’ and these exact
solutions could tolerate some errors. The Jacobi iteration method is an easy
way to obtain approximate solutions [8, 9], and it mainly involves matrix-
vector multiplications and vector additions, which are of fine-grained
parallelism. An alternative method is to use sparse approximate inverses
(SAI), which is based on the decay of inverses of sparse matrices [11, 12,
19]. With SAIs, the two triangular solvers in ILU preconditioning can be
replaced by two matrix-vector multiplications. There are several ways to
compute the SAI of a matrix, for example Frobenius norm minimization and
incomplete inverse factorization [5, 6, 7, 20]. Also, there are some SAIs
specially designed for triangular systems [4, 16, 22]. In [13, 21], the
truncated Neumann expansions play a similar role to the SAIs of ILU factors.
Some recent processes of this topic are introduced in [1].

When deriving the method we propose in this paper, we compress the Jacobi
iteration method to a series firstly. As the iterative matrix in Jacobi
iteration method for triangular systems is strictly triangular, the series is
finite. We take the truncation of this series as the approximate inverse of
triangular matrix and use some dropping strategies to make the inverse keep
enough sparsity. We rewrite the truncation of the series as a recursive
formula, by which we can avoid computing and storing high-order matrix powers.
Then the SAIT algorithm is written in a quite brief way and it mainly involves
sparse matrix-matrix multiplications and some dropping rules. The SAIT
matrices can be applied to ILU preconditioning. The two triangular solvers can
be replaced by two matrix-vector multiplications with a pair of SAIT matrices,
which can be fine-grained parallelized.

This paper is organized as follows. In section 2, we propose an exact inverse
for triangular matrices base on Jacobi iteration and then construct the
general formulation of the SAIT algorithm. In section 3, we use two dropping
strategies in the SAIT algorithm. In section 4, we take the SAIT matrices as
preconditioners in solving linear systems and eigenvalue problems with
preconditioned iterative methods. In section 5, there are some conclusions.

##  2 The basic idea of the SAIT algorithm

Let $T\in\mathbb{R}^{n\times n}$ be a triangular matrix. We solve the linear
equation

$\begin{split}Tx=b\end{split}$

by Jacobi iteration method

$\begin{split}x^{k}=D^{-1}(D-T)x^{k-1}+D^{-1}b.\end{split}$ (1)

Here, $D$ is the diagonal matrix of $T$. For the sake of simplicity, we denote
${T}_{0}=I-D^{-1}T$. ${T}_{0}$ is a strictly triangular matrix. Then, (1)
becomes

$\begin{split}x^{k}={T}_{0}x^{k-1}+D^{-1}b.\end{split}$

Letting the initial datum $x^{0}=0$, we list the solution of each iteration:

$\begin{split}x^{1}&=D^{-1}b\\\
x^{2}&={T}_{0}x^{1}+D^{-1}b=({T}_{0}+I)D^{-1}b\\\
x^{3}&={T}_{0}x^{2}+D^{-1}b=({T}_{0}^{2}+{T}_{0}+I)D^{-1}b\\\ &\cdots\\\
x^{k}&={T}_{0}x^{k-1}+D^{-1}b=\left(\sum_{i=0}^{k-1}{T}_{0}^{i}D^{-1}\right)b.\\\
\end{split}$

Denoting

$\begin{split}M_{k}=\sum_{i=0}^{k-1}{T}_{0}^{i}D^{-1},\end{split}$ (2)

then we have

$\begin{split}x^{k}=M_{k}b.\end{split}$

As ${T}_{0}$ is a strictly triangular matrix, we have ${T}_{0}^{k}=0$ when
$k\geq n$. Then we know that the exact inverse of $T$ is a finite series of
(2), i.e.

$\begin{split}T^{-1}=M_{n}\equiv\sum_{i=0}^{n-1}{T}_{0}^{i}D^{-1}.\end{split}$
(3)

In the following, we construct approximate inverses for the triangular matrix
$T$ based on its exact inverse (3). The approximation comes from two aspects,
truncation and dropping. Let we study a truncated term $M_{m}$ of (3), where
$m<n-1$ and

$\begin{split}M_{m}=D^{-1}+{T}_{0}D^{-1}+{T}_{0}^{2}D^{-1}+\cdots+{T}_{0}^{m}D^{-1}.\end{split}$
(4)

In order to avoid computing and storing the high-order terms ${T}_{0}^{i}$
$(i=1,\cdots,m)$ in (4), we rewrite $M_{k}$ as an equivalent form:

$\begin{split}M_{0}D&=I,\\\ M_{1}D&={T}_{0}+I,\\\
M_{2}D&={T}_{0}({T}_{0}+I)+I,\\\ M_{3}D&={T}_{0}({T}_{0}({T}_{0}+I)+I)+I,\\\
&\cdots\\\ \end{split}$ (5)

This is the Horner’s method in calculating the values of polynomial functions.
The same method was also proposed by a Chinese mathematician _Qin Jiushao_ in
$13^{th}$ century Song dynasty [23]. Furthermore, we present (5) in a
recursive formula:

$\begin{split}M_{0}D&=I,\\\ M_{1}D&={T}_{0}(M_{0}D)+I,\\\
M_{2}D&={T}_{0}(M_{1}D)+I,\\\ M_{3}D&={T}_{0}(M_{2}D)+I,\\\ &\cdots\\\
M_{m}D&={T}_{0}(M_{m-1}D)+I.\\\ \end{split}$ (6)

To make the matrix more sparse, we perform dropping rules after each
multiplication by ${T}_{0}$ in (6). This strategy is summarized in Algorithm
1.

Algorithm 1 Sparse Approximate Inverse for Triangular matrices base on Jacobi
iteration (SAIT)

1: Set the initial data $M=I$, let $D$ be the diagonal matrix of $T$ and let
${T}_{0}=I-D^{-1}T$.

2: for $k=1,2,\cdots,m$ do

3: $M={T}_{0}M+I$

4: execute dropping rules on $M$

5: end for

6: $M=MD^{-1}$

##  3 Dropping strategies

We use two dropping strategies, threshold-based and pattern-based, to
concretize Algorithm 1.

Threshold-based. In step 3 of Algorithm 1, we drop the entries whose
magnitudes are small than $\tau$ $(\tau<1)$ in $M$. As the diagonal entries of
$M$ are all units, the setting $(\tau<1)$ can guarantee the matrix $M$ being
always full-rank during the iterations. This strategy is presented in
Algorithm 2. We paste the Matlab code for this algorithm in Appendix A.

Algorithm 2 Threshold-based Sparse Approximate Inverse for Triangular matrices
SAIT_Thr$(\tau,m)$

1: Set the initial data $M=I$, let $D$ be the diagonal matrix of $T$ and let
${T}_{0}=I-D^{-1}T$.

2: for $k=1,2,\cdots,m$ do

3: $M={T}_{0}M+I$

4: drop the entries in $M$ whose magnitudes are small than $\tau$ $(\tau<1)$

5: end for

6: $M=MD^{-1}$

Pattern-based. We set a fixed sparsity pattern ${\boldsymbol{S}}$ in advance,
and drop the entries out of ${\boldsymbol{S}}$ in each iteration. The patterns
of the power functions of the original matrix are often used in computing
sparse approximate inverses [20]. We study the sparsity pattern of $T^{p}$ by

$\begin{split}T^{p}=((T-D)+D)^{p}=\sum_{i=0}^{p}C_{p}^{i}D^{p-i}(T-D)^{i}.\end{split}$

Since $(T-D)$ shares the same pattern with ${T}_{0}$, compared with (2), we
know that $T^{p}$ has the same pattern with $M_{p+1}$. As $M_{p}$ tends to the
exact inverse with $p\to n$, this pattern can be regarded as a truncated
pattern of the exact inverse. When $p=0$, it is the Jacobi preconditioner,
while when $p=1$, it is the pattern of the original matrix $T$. The pattern of
$T^{p}$ can be obtained by Algorithm 1 without dropping any entries in the
first $p$ iterations. Then, we keep the pattern of $T^{p}$, and in the
following iterations, drop the entries out of this pattern. The number $p$
should not be too large, as the number of nonzeros in the pattern of $M_{p}$
increases rapidly with $p$ increasing. This dropping strategy is summarized in
Algorithm 3.

Algorithm 3 Pattern-based Sparse Approximate Inverse for Triangular matrices
SAIT_Pat$(p,m)$

1: Set the initial data $M=I$, let $D$ be the diagonal matrix of $T$ and let
${T}_{0}=I-D^{-1}T$.

2: for $k=1,2,\cdots,p$ do

3: $M={T}_{0}M+I$

4: end for

5: keep the sparsity pattern ${\boldsymbol{S}}$ of $M$

6: for $k=1,2,\cdots,m$ do

7: $M={T}_{0}M+I$

8: drop the entries in $M$ out of the pattern ${\boldsymbol{S}}$

9: end for

10: $M=MD^{-1}$

###  3.1 Balance between the accuracy and iteration counts

Let $L$ and $U$ be the two triangular factors of incomplete LU (ILU)
factorization of a matrix $A$, and let $M_{L}$ and $M_{U}$ be the SAIT
matrices for them respectively. In ILU preconditioning, the multiplication
$LU$ can be regarded as an approximation of $A$, i.e. $LU\approx A$. If we
take $M_{L}$ and $M_{U}$ as the approximations of $L^{-1}$ and $U^{-1}$,
respectively, the multiplication $M_{U}M_{L}$ can be take an approximate
inverse of $A$, i.e.

$\begin{split}M_{U}M_{L}\approx U^{-1}L^{-1}\approx A^{-1}.\end{split}$

When the ILU is taken as a preconditioner in an iterative method, we replace
the two triangular equations $U(Lx)=y$ by two matrix-vector multiplications
$x=M_{U}(M_{L}y)$.

In the Algorithm 2 and 3, if we use smaller threshold $\tau$ and larger $p$,
we can obtain more accurate approximate inverses, since the SAIT matrices
contains more nonzero entries. However, higher accuracy of SAIT matrices do
not mean shorter runtime in solving a system with this preconditioner.

We list the following terms related to an iterative solver with the SAIT
preconditioners:

  * •

$iter(M_{L},M_{U})$: the iteration count of an iterative solver with $M_{L}$
and $M_{U}$ in preconditioning procedure.

  * •

$comp(M_{L},M_{U})$: the runtime of computing the two matrix-vector
multiplications $M_{U}(M_{L}y)$.

  * •

$othercomp$: the runtime of computing other terms in each iteration except the
preconditioning procedure.

Then, the runtime of a preconditioned iterative method can be estimated
roughly by the following formula:

$\begin{split}runtime=iter(M_{L},M_{U})\times\bigg{(}comp(M_{L},M_{U})+othercomp\bigg{)}.\end{split}$
(7)

We define the ratios of a SAIT matrix $M_{T}$ compared with its original
matrix $T$:

$\begin{split}r=\frac{nnz(M_{T})}{nnz(T)},\end{split}$

where $nnz(\cdot)$ denotes the number of nonzeros in a matrix. This value is
to illustrate the number of nonzeros in a SAIT matrix, and it correlates to
the runtime of the part $comp(M_{L},M_{U})$.

When we use iterative methods to solving a linear system, if the iterative
method and the computer environment are fixed, the term $othercomp$ is usually
fixed. If we use more accurate SAIT matrices, the iteration count
$iter(M_{L},M_{U})$ decreases. However, we can not expect that the SAIT
preconditioners can reduce the iteration count less than exact triangular
solvers. It is probable that after exceeding some point, more nonzeros in
$M_{L}$ and $M_{U}$ can not reduce the $iter(M_{L},M_{U})$ further. On the
other side, more nonzeros can cause the term $comp(M_{L},M_{U})$ increasing,
and sequentially, the total $runtime$ increasing.

In order to reduce the runtime, it should make a balance between the accuracy
of SAIT matrices and the iteration count. As the accuracy is mainly decided by
the threshold $\tau$ and the number $p$ in the Algorithm 2 and 3,
respectively, there should be some $\tau$ or $p$ that make the runtime of an
iterative solver with SAIT shortest. We call such parameters _Optimal
parameters_. The optimal parameters may vary with different problems,
different methods and different computer environments. We will verify the
discussion above using numerical experiments in section 4.3.

##  4 Numerical experiments

In the numerical experiments, the main test models are the three-dimensional
Laplace equation and its eigenvalue problem with homogeneous Dirichlet
boundary condition:

$\begin{split}-\Delta u&=f\quad\,\,\;\text{in}\,\,\Omega,\qquad
u=0\quad\text{on}\,\,\partial\Omega\\\ \text{and}\quad-\Delta u&=\lambda
u\quad\text{in}\,\,\Omega,\qquad
u=0\quad\text{on}\,\,\partial\Omega,\end{split}$ (8)

where $\Omega=[0,1]^{3}$. We use a $102\times 102\times 102$ uniform mesh and
finite difference method to discretize them. We obtain the corresponding
matrix problems

$\begin{split}Au_{h}&=f_{h}\\\ \text{and}\quad
Au_{h}&=\lambda_{h}u_{h}.\end{split}$ (9)

Here, $A$ is a square matrix with $10^{6}$ rows/columns and it is symmetric
positive definite (SPD). We use preconditioned conjugate gradient method (PCG)
and LOBPCG method [17] to compute the linear equations and eigenvalue problem
in (9), respectively.

We take the level-0 and level-1 ILU factors of the matrix $A$ as the original
triangular systems and use the SAIT algorithms to approach their inverses.
Then, we take the SAIT matrices as the preconditioners in PCG and LOBPCG. The
SAIT matrices are generated by different SAIT algorithms with different
parameters. We test and compare their effects in preconditioning the iterative
solvers.

The code is implemented in Matlab and is run on a laptop with an Intel
i7-6700HQ CPU with 16 GB RAM and a GTX970M GPU with 3 GB memory. We take
$10^{-10}$ as the uniform stopping criteria for all the numerical experiments.

###  4.1 SAIT_Thr$(\tau,m)$ and SAIT_Pat$(p,m)$

We test the numbers of nonzeros in the SAIT matrices generated by the two
dropping strategies, the SAIT_Thr$(\tau,m)$ in Algorithm 2 and SAIT_Pat$(p,m)$
in Algorithm 3. Figure 1 shows the ratios of nonzeros in SAIT matrices
generated by SAIT_Thr$(\tau,m)$ with different threshold
$\tau=0,\;0.01,\;0.02,\;0.05$ and increasing iteration number $m$. The setting
$\tau=0$ means that there is no entry dropped in each iteration. In this case,
the numbers of nonzeros increase rapidly. For other nonzero $\tau$, the
numbers of nonzeros reach different stable states after several iterations.
The smaller $\tau$ results in the more nonzeros in its corresponding SAIT
matrix. For SAIT_Pat$(p,m)$, a fixed $p$ leads to a fixed pattern, which means
that the number of nonzeros of its SAIT matrix is fixed. We put the ratios of
the SAIT matrices generated by this algorithm in the legends of Figure 3.

Figure 1: The ratios of the nonzeros in the SAIT matrices generated by
SAIT_Thr$(\tau,m)$ with different $\tau$ and $m$ for level-0 (left) and
level-1 (right) ILU factors.

Next, we test the effects of SAIT matrices as preconditioners in PCG method.
For the matrices generated by SAIT_Thr$(\tau,m)$, the iteration counts of PCG
decrease with increasing $m$ until reaching different stable states, shown in
Figure 2. When using SAIT_Pat$(p,m)$, the iteration counts keep stable with
$m$ varying, shown in Figure 3. For the both algorithms, the small threshold
$\tau$ or larger $p$ (more nonzeros) means less iteration count. However,
there is no case that the SAIT preconditioners can reduced the iteration
counts less than the exact triangular solvers.

Figure 2: The iteration counts of PCG with SAIT preconditioners generated by
SAIT_Thr$(\tau,m)$ with different $\tau$ and $m$ for level-0 (left) and
level-1 (right) ILU factors.

Figure 3: The iteration counts of PCG with SAIT preconditioners generated by
SAIT_Pat$(p,m)$ with different $p$ and $m$ for level-0 (left) and level-1
(right) ILU factors.

We also compare the two dropping strategies with a fixed parameter $m=10$. As
it is shown in Figure 4, with the same numbers of nonzeros, the PCG iteration
counts of SAIT_Thr$(\tau,10)$ are less than that of SAIT_Pat$(p,10)$.

Figure 4: The comparisons of the iteration counts of PCG with SAIT
preconditioner generated by SAIT_Thr$(\tau,m)$ and SAIT_Pat$(p,m)$ for level-0
(left) and level-1 (right) ILU factors.

###  4.2 The balance

Table 1 shows the iteration counts and runtime in solving the linear equation
in (9) using PCG with the SAIT preconditioners. We observe that less iteration
count does not means shorter runtime. We test more threshold parameters in
SAIT_Thr$(\tau,m)$, shown in Figure 5. The iteration count does not decrease
continuously with smaller threshold $\tau$ (more nonzeros) after some point.
Even though we use the exact inverses which are usually full triangular
matrices, the iteration count is equal to the exact triangular solver, which
is not an arbitrarily small number. However, a SAIT matrix with too many
nonzeros can increase the cost in computing matrix-vector multiplications.
From the relation between the runtime and numbers of nonzeros, shown in the
right picture of Figure 5, there should exists an optimal threshold $\tau$ in
respect of the runtime. Such optimal parameter can be found by numerical
experiments for each specific problem.

Figure 5: Left: the relation between the PCG iteration counts and the numbers
of nonzeros in the SAIT preconditioners. Right:the relation between the
runtime of PCG method and the numbers of nonzeros in the SAIT preconditioners.

###  4.3 Comparison with Jacobi iteration

Since the SAIT algorithms are designed based on Jacobi iteration method, if we
use the Jacobi method in the preconditioning procedure directly, it is
equivalent to the SAIT matrices generated by SAIT_Thr$(0,m)$, which are more
accurate than the cases with thresholds $\tau>0$. In addition, it requires
less memory, since a SAIT matrix usually has more nonzeros than its original
matrix. We list the iteration counts and runtime of the Jacobi method in Table
2. The PCG iteration counts with Jacobi method can be almost the same as the
exact solver after several iteration. However, its runtime is much longer than
that of SAIT preconditioners in Table 1. The reason is that Jacobi method
involves several matrix-vector multiplications and vector-vector additions in
each PCG iteration, while there are only two matrix-vector multiplications
with the SAIT preconditioners.

| level-0 ILU | level-1 ILU
---|---|---
ratio | iter | time (s) | ratio | iter | time (s)
SAIT_Thr(0.05,10) | 1.74 | 189 | _9.95_ | 1.00 | 184 | _9.97_
SAIT_Thr(0.02,10) | 2.73 | 168 | 10.52 | 3.37 | 133 | 12.73
SAIT_Thr(0.01,10) | 4.92 | _154_ | 13.24 | 5.17 | _121_ | 15.48
SAIT_Pat(1,10) | 1.00 | 228 | 10.83 | 1.00 | 229 | 12.22
SAIT_Pat(2,10) | 2.48 | 177 | 10.95 | 3.25 | 158 | 15.37
SAIT_Pat(3,10) | 4.92 | _154_ | 13.25 | 7.54 | 129 | 22.09
Table 1: The ratios the nonzeros, iteration count and runtime in solving the 3D FDM system using PCG method with different SAIT preconditoners. The red italic numbers are the least iteration account or the shortest runtime of the SAIT preconditoners generated by level-0 and level-1 ILU factors, respectively. Jacobi iter | level-0 ILU | level-1 ILU
---|---|---
iter | time (s) | iter | time (s)
1 | 423 | 29.88 | 423 | 31.87
2 | 229 | 21.03 | 240 | 24.78
3 | 173 | _19.33_ | 169 | 22.10
4 | 152 | 20.43 | 134 | _21.19_
5 | 152 | 23.23 | 116 | 21.48
6 | 149 | 25.94 | 108 | 23.45
7 | 147 | 28.19 | 104 | 24.74
8 | 146 | 30.76 | 101 | 27.21
9 | 145 | 33.67 | 97 | 28.38
10 | 145 | 36.60 | 98 | 31.31
11 | 145 | 39.53 | 98 | 33.90
12 | 145 | 42.28 | 98 | 36.50
13 | 145 | 45.66 | 98 | 39.25
14 | 145 | 48.44 | 97 | 41.43
15 | 145 | 51.54 | 97 | 47.38
Table 2: The iteration counts and runtime in solving the 3D FDM system in (9)
using PCG method with Jacobi iteration method as preconditioner. The red
italic numbers are the shortest runtime of level-0 and level-1 ILU factors,
respectively.

###  4.4 Acceleration by GPU

The main computations of the basic Conjecture Gradient method are matrix-
vector multiplications and vector-inner products. With the SAIT
preconditioners, the preconditioning procedure are two matrix-vector
multiplications. These operations can be highly parallelized. We use the
function $gpuArray$ in Matlab to input the data into GPU. We run the solver
with the same preconditioners in Table 1 on the GPU. Though the iteration
counts are the similar with the experiments on CPU, the codes are accelerated
several times by GPU, shown in Table 3.

| level-0 ILU | level-1 ILU
---|---|---
ratio | iter | time (s) | ratio | iter | time (s)
SAIT_Thr(0.05,10) | 1.74 | 189 | 1.34 | 1.00 | 184 | 1.34
SAIT_Thr(0.02,10) | 2.73 | 168 | 1.33 | 3.37 | 133 | 1.59
SAIT_Thr(0.01,10) | 4.92 | 154 | 1.81 | 5.17 | 121 | 1.81
SAIT_Pat(1,10) | 1.00 | 228 | 1.47 | 1.00 | 229 | 1.67
SAIT_Pat(2,10) | 2.48 | 177 | 1.36 | 3.25 | 158 | 1.92
SAIT_Pat(3,10) | 4.92 | 154 | 1.80 | 7.54 | 129 | 2.38
Table 3: The same codes in Table 1 are run on GPU.

###  4.5 Preconditioned solver for eigenvalue problems

When we use the simultaneous preconditioned method LOBPCG to compute the
eigenvalue problem in (9), there are many vectors to be dealt with in the
preconditioning procedure of each iteration. With SAIT preconditioners, this
procedure can be done through two matrix-matrix multiplications. Here, we use
the LOBPCG method to compute the first 4 eigenvalues in this problems. We use
level-0 and level-1 ILU factorizations, and use SAIT_Thr$(0.01,10)$ to
generated the SAIT matrices. Figure 6 and 7 show their results, respectively.

Figure 6: The convergence history in solving Laplace eigenvalue problem using
LOPBCG method with the exact solver and a SAT preconditioner generate by
SAIT_Thr$(0.03,10)$ for level-0 ILU factors in the preconditioning procedure.

Figure 7: The convergence history in solving Laplace eigenvalue problem using
LOPBCG method with the exact solver and a SAT preconditioner generated by
SAIT_Thr$(0.03,10)$ for level-1 ILU factors in the preconditioning procedure.

###  4.6 More examples

We apply the SAIT algorithms to more examples, shown in Table 4. These
examples are from the University of Florida sparse matrix collection [10].
Table 5 are the results using the SAIT_Thr$(\tau,10)$ preconditioners. Here,
the ratios of nonzeros in the SAIT matrices generated by their corresponding
thresholds (may be not optimal) are mainly between 1 and 2. We also try
SAIT_Pat$(p,10)$ for these examples, shown in Table 6. Comparing with the
results in the two tables, we find that the threshold-based
SAIT_Thr$(\tau,10)$ preconditioners preform better than the pattern-based
SAIT_Pat$(p,10)$ for these problems. In Table 5, we also compare the iteration
counts of SAIT preconditioners with the exact solver. Even though more
iterations, the runtime of SAIT can be reduced by parallelization as it is
discussed in section 4.4.

| row/column |  nnz
---|---|---
apache1 | 80800 | 542184
apache2 | 715176 | 4817870
thermal1 | 82654 | 574458
thermal2 | 1228045 | 8580313
parabolic_fem | 525825 | 3674625
G3_circuit | 1585478 | 7660826
ecology2 | 999999 | 4995991
thermomech_dM | 204316 | 1423116
thermomech_TC | 102158 | 711558
Table 4: The numbers of rows and nonzeros of some SPD matrices in [10]. | no p.c. | level-0 ILU | level-1 ILU
---|---|---|---
| exact | $\tau$ | ratio | SAIT |  | exact | $\tau$ | ratio | SAIT |
apache1 | 3777 | 365 | 0.05 | 1.46 | 439 | 20.3% | 249 | 0.03 | 1.80 | 316 | 26.9%
apache2 | 5528 | 882 | 0.05 | 1.39 | 1092 | 23.8% | 587 | 0.03 | 1.70 | 797 | 35.8%
thermal1 | 1707 | 651 | 0.05 | 1.43 | 703 | 8.0% | 363 | 0.04 | 1.76 | 435 | 19.8%
thermal2 | 6626 | 2555 | 0.05 | 1.42 | 2763 | 8.1% | 1401 | 0.04 | 1.72 | 1674 | 19.5%
parabolic_fem | 3515 | 1423 | 0.1 | 0.96 | 1640 | 15.2% | 845 | 0.04 | 1.32 | 946 | 12.0%
thermomech_dM | 89 | 10 | 0.06 | 1.01 | 11 | 10.0% | 6 | 0.02 | 1.02 | 8 | 33.3%
thermomech_TC | 89 | 10 | 0.06 | 1.01 | 11 | 10.0% | 6 | 0.02 | 1.02 | 8 | 33.3%
ecology2 | 7127 | 2123 | 0.08 | 2.00 | 2830 | 33.3% | 1303 | 0.06 | 3.25 | 1942 | 49.0%
G3_circuit | 21205 | 1182 | 0.08 | 2.07 | 1582 | 33.8% | 643 | 0.08 | 2.38 | 1174 | 82.6%
Table 5: The iteration counts of PCG method with the exact solver and SAIT preconditioners generated by SAIT_Thr$(\tau,10)$. | level-0 ILU | level-1 ILU
---|---|---
| m = 2 | m = 3 | m = 2 | m = 3
| ratio | SAIT | ratio | SAIT | ratio | SAIT | ratio | SAIT
thermomech_dM | 1 | 11 | 1.60 | 10 | 1 | 10 | 2.37 | 6
thermomech_TC | 1 | 11 | 1.60 | 10 | 1 | 10 | 2.37 | 6
thermal1 | 1 | 847 | 1.85 | 690 | 1 | 743 | 2.49 | 440
thermal2 | 1 | 3284 | 1.85 | 2674 | 1 | 2868 | 2.48 | 1635
parabolic_fem | 1 | 1678 | 1.56 | 1451 | 1 | 1444 | 2.38 | 893
apache1 | 1 | 3252 | 2.41 | 2236 | 1 | 4102 | 3.11 | 2570
apache2 | 1 | 2753 | 2.42 | 1818 | 1 | 3521 | 3.12 | 2500
G3_circuit | 1 | 3993 | 2.04 | 2751 | 1 | 3713 | 2.37 | 2245
ecology2 | 1 | 3738 | 2.00 | 2799 | 1 | 3765 | 2.25 | 2549
Table 6: The iteration counts of PCG method with the SAIT preconditioners
generated by SAIT_Pat$(p,10)$.

##  5 Conclusion

We derive an exact inverse for triangular matrix trough Jacobi iteration
method. The inverse is a finite series. We take the truncation of this series
as the approximate inverse for a triangular matrix. To make the approximate
inverses more sparse, we propose two dropping strategies. We apply the SAIT
matrices to iterative method with ILU preconditioners. Then the two triangular
solvers in ILU preconditioning are replaced by two matrix-vector
multiplications, which can be fine-grained parallelized. We present some
numerical examples to show the effects of SAIT preconditioners. Though the
iteration counts of SAIT matrices are larger than the exact triangular
solvers, the runtime can be reduced by parallelization. In this paper, we only
present the basic study on this method. There are many things about this
method are remained to be studied further, for example, more efficient
dropping strategies, choices of the parameters and parallel implementations.

##  Appendix A The Matlab code of SAIT_Thr$(\tau,m)$

⬇

function M = SAIT_Thr(T,tau,m)

I = speye(size(T,1));

Dn = diag(diag(T).^(-1));

T0 = I - Dn*T;

M = I;

for i = 1 : m

M = T0*M + I;

P = (abs(M)>tau);

M = M.*P;

end

M = M*Dn;

end

## References

  * [1] Ahmad Abdelfattah, Hartwig Anzt, Jack Dongarra, Mark Gates, Azzam Haidar, Jakub Kurzak, Piotr Luszczek, Stanimire Tomov, Ichitaro Yamazaki, and Asim YarKhan.  Linear algebra software for large-scale accelerated multicore computing.  Acta Numerica, 25:1–160, 2016.
  * [2] Edward Anderson and Youcef Saad.  Solving sparse triangular linear systems on parallel computers.  Internat. J. High Speed Comput., 1(01):73–95, 1989.
  * [3] Hartwig Anzt, Edmond Chow, and Jack Dongarra.  ParILUT—a new parallel threshold ILU factorization.  SIAM J. Sci. Comput., 40(4):C503–C519, 2018.
  * [4] Hartwig Anzt, Thomas K. Huckle, Jürgen Bräckle, and Jack Dongarra.  Incomplete sparse approximate inverses for parallel preconditioning.  Parallel Comput., 71:1–22, 2018.
  * [5] Michele Benzi.  Preconditioning techniques for large linear systems: a survey.  J. Comput. Phys., 182(2):418–477, 2002.
  * [6] Michele Benzi, Carl D. Meyer, and Miroslav Tůma.  A sparse approximate inverse preconditioner for the conjugate gradient method.  SIAM J. Sci. Comput., 17(5):1135–1149, 1996.
  * [7] Michele Benzi and Miroslav Tuma.  A comparative study of sparse approximate inverse preconditioners.  Appl. Numer. Math., 30(2-3):305–340, 1999.
  * [8] Edmond Chow, Hartwig Anzt, Jennifer Scott, and Jack Dongarra.  Using Jacobi iterations and blocking for solving sparse triangular systems in incomplete factorization preconditioning.  J. Parallel Distrib. Comput., 119:219–230, 2018.
  * [9] Edmond Chow and Aftab Patel.  Fine-grained parallel incomplete LU factorization.  SIAM J. Sci. Comput., 37(2):C169–C193, 2015.
  * [10] Timothy A. Davis and Yifan Hu.  The University of Florida sparse matrix collection.  ACM T. Math. Software (TOMS), 38(1):1, 2011.
  * [11] Stephen Demko, William F. Moss, and Philip W Smith.  Decay rates for inverses of band matrices.  Math. Comp., 43(168):491–499, 1984.
  * [12] Victor Eijkhout and Ben Polman.  Decay rates of inverses of banded M-matrices that are near to toeplitz matrices.  Linear Algebra Appl., 109:247–277, 1988.
  * [13] Ivar Gustafsson and Gunhild Lindskog.  Completely parallelizable preconditioning methods.  Numer. Linear Algebra Appl., 2(5):447–465, 1995.
  * [14] Pascal Hénon and Yousef Saad.  A parallel multistage ILU factorization based on a hierarchical graph decomposition.  SIAM J. Sci. Comput., 28(6):2266–2293, 2006.
  * [15] David Hysom and Alex Pothen.  A scalable parallel algorithm for incomplete factor preconditioning.  SIAM J. Sci. Comput., 22(6):2194–2215, 2001.
  * [16] Carlo Janna, Massimilano Ferronato, and Giuseppe Gambolati.  A block FSAI-ILU parallel preconditioner for symmetric positive definite linear systems.  SIAM J. Sci. Comput., 32(5):2468–2484, 2010.
  * [17] Andrew V. Knyazev.  Toward the optimal preconditioned eigensolver: Locally optimal block preconditioned conjugate gradient method.  SIAM Journal on Scientific Computing, 23(2):517–541, 2001.
  * [18] Mardochée Magolu monga Made and Henk A. van der Vorst.  A generalized domain decomposition paradigm for parallel incomplete LU factorization preconditionings.  Future Gener. Comp. Sy., 17(8):925–932, 2001.
  * [19] Reinhard Nabben.  Decay rates of the inverse of nonsymmetric tridiagonal and band matrices.  SIAM J. Matrix Anal. Appl., 20(3):820–837, 1999.
  * [20] Yousef Saad.  Iterative methods for sparse linear systems, volume 82.  SIAM, 2003.
  * [21] Henk A. van der Vorst.  A vectorizable variant of some ICCG methods.  SIAM J. Sci. Statist. Comput., 3(3):350–356, 1982.
  * [22] Arno C. N. van Duin.  Scalable parallel preconditioning with the sparse approximate inverse of triangular matrices.  SIAM J. Matrix Anal. Appl., 20(4):987–1006, 1999.
  * [23] Wen-tsun Wu.  Grand Series of Chinese Mathematics (in Chinese), volume V.  Beijing Normal University Publishing House, 2000.