# On the origin of observed cosmic ray spectrum below 100 TV

Mikhail A. Malkov<EMAIL_ADDRESS>Department of Physics and CASS, University
of California San Diego, La Jolla, CA 92093, USA  Igor V. Moskalenko
<EMAIL_ADDRESS>Hansen Experimental Physics Laboratory and Kavli Institute
for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA
94305, USA

###### Abstract

Significant progress in cosmic ray (CR) studies was achieved over the past
decade. Particularly important are precise measurements of primary and
secondary species in the TV rigidity domain that show a bump in the spectra of
CR species from 0.5–50 TV. In this letter, we argue that it is likely caused
by a stellar bow- or wind-termination shock that reaccelerates preexisting
CRs, which further propagate to the Sun along the magnetic field lines. This
single universal process is responsible for the observed spectra of all CR
species in the rigidity range below 100 TV. A viable candidate is Epsilon
Eridani star at 3.2 pc from the Sun, which is well-aligned with the direction
of the local magnetic field. We provide a simple formula that reproduces the
spectra of all CR species with only two nonadjustable shock parameters,
uniquely derived from the proton data. We show how our formalism predicts
helium and carbon spectra and the B/C ratio.

## Discovery of the CR “Bump.”—

The last decade was marked by discoveries of new features in the spectra of CR
species. Unexpected breaks and variations in spectral indices are being
unveiled with high confidence in the GV–TV rigidity range that deemed as well-
studied [ATIC06, 2010ApJ...714L..89A, 2011Sci...332...69A,
2009BRASP..73..564P, 2011ApJ...728..122Y, 2014PhRvL.112o1103A,
2015PhRvL.114q1103A, 2015PhRvL.115u1101A, 2017PhRvL.119y1101A,
2018PhRvL.120b1101A, 2018PhRvL.121e1103A, PhysRevLett.124.211102,
2018JETPL.108....5A, 2019ARep...63...66A, 2019AdSpR..64.2546G,
2019PhRvL.122r1102A, eaax3793, Aguilar2021, 2019AdSpR..64.2559G]). These
features bear the signatures of CR acceleration processes and their
propagation history.

First discovered was a flattening above $R_{{\rm br}}^{\prime}\approx 0.5$ TV
in rigidity in the CR proton and helium spectra [ATIC06, 2010ApJ...714L..89A,
2011Sci...332...69A], which propelled a slew of interpretations. However, with
improved measurements of spectra of several primary ($p$, He, C, O, Ne, Mg,
Si) and secondary111The mostly secondary species are those, which are rare in
CR sources, but abundant in CRs. They are produced primarily in fragmentations
of heavier species over their lifetime in the Galaxy (Li, Be, B) species by
AMS-02 [2015PhRvL.114q1103A, 2015PhRvL.115u1101A, 2017PhRvL.119y1101A,
2018PhRvL.120b1101A, 2018PhRvL.121e1103A, PhysRevLett.124.211102], only two
interpretations remain feasible. These are (i) the intrinsic spectral break in
the CR injection spectrum or populations of sources with soft and hard
spectra, a so-called “injection scenario” [2012ApJ...752...68V], and (ii) a
break in the spectrum of interstellar turbulence translated into the break in
the diffusion coefficient, a so-called “propagation scenario”
[2012ApJ...752...68V, 2012PhRvL.109f1101B, AloisioKink_2015]. The latter has
very few free parameters and is well-supported [2020ApJ...889..167B,
Boschini_2020] by the AMS-02 data.

A turning point was the recent discovery of a steepening in the proton and
helium spectra at $R_{{\rm br}}^{\prime\prime}\sim$ 13 TV [ATIC06,
2018JETPL.108....5A, 2019ARep...63...66A, 2019AdSpR..64.2546G,
2019PhRvL.122r1102A, eaax3793] (see parameterizations in [Boschini_2020]).
Together with the lower rigidity break $R_{{\rm br}}^{\prime}$, the two
features likely result from a single physical process that generates a “bump,”
in the CR spectrum between 0.5–50 TV. It peaks by a factor of 2.5–3 near
$\approx$13 TV above the background CR spectrum.

The discovery of the second break challenges conventional interpretations of
the CR spectra outlined above. This finding also refutes earlier hypotheses of
a contribution from a nearby CR accelerator, such as a supernova remnant
(SNR), that accelerates the bump particles out of the interstellar medium
(ISM) [e.g., 2020ApJ...903...69F, Fornieri2020, 2020FrPhy..1624501Y], because
such scenario implies no breaks in the secondary component
[2012ApJ...752...68V]. Moreover, the unprecedented accuracy with which the
sharp break at $R_{{\rm br}}^{\prime}$ in the spectrum is measured severely
constrains its remote or global origins as they would lead to a smoother
transition.

A clue to the origin of the bump is the break in the spectra of secondary
species that are also flattening above the first break ($R_{{\rm
br}}^{\prime}$), similarly to the primaries, albeit with a different spectral
index. It implies that the secondary species have already been present in the
CR mixture that the bump is made of. This leads us to conclude that (i) the
bump has to be made out of the preexisting CRs with all their primaries and
secondaries that have spent millions of years in the Galaxy. Meanwhile, the
sharpness of the breaks indicates that (ii) the bump is formed locally, a
short time before we observed it.

## Local CR Reacceleration.—

The consequences of the above reasoning are two-fold. First, the bump must
consist of reaccelerated background CRs. A salient feature of diffusive
reacceleration in weak shocks is that stronger reacceleration experience CR
elements with steeper spectra. Second, a search for a local source led us to
an idea to look for a bow- or wind-termination shock of a passing star nearby.
Estimates show that such a shock should be within 3–10 pc of the Sun
[MalkMosk_2021]. In considering the candidates, one encounters an intriguing
alignment of the local magnetic field with the Epsilon Eridani star at 3.2 pc
from the Sun, which has an extensive 8000 au astrosphere. Nevertheless, other
stars/shocks cannot be ruled out by the current data.

## Model Scrutiny.—

We employ a nearby shock reacceleration model, complemented with CR
propagation along a magnetic flux tube pointing to the Solar System. It
depends on an unknown shock Mach number and its size-distance relation
[MalkMosk_2021]. The key advantage of the model is that these parameters are
uniquely derived from the best measured CR proton spectrum or its substitute.
As the latter, in [MalkMosk_2021] we used “synthetic” data provided in
[Boschini_2020], i.e. a physically motivated local interstellar spectrum (LIS)
of CR protons tuned to the data of different instruments with adjustments for
the appropriate heliospheric modulation.

As the slope and normalization of CR background LIS [Boschini_2020] below the
first break are now accurately derived, the model takes them as fixed input
parameters. It then fully predicts the entire bump structure including two
break rigidities ($R_{\rm br}^{\prime}$, $R_{\rm br}^{\prime\prime}$), two
respective jumps in the spectral index and two curvatures of the spectrum
(sharpness) at the breaks (altogether six parameters).

The proton spectrum fit is critically determined by the turbulence index $s$
found to be 3/2 (Iroshnikov-Kraichnan spectrum [1964SvA.....7..566I,
1965PhFl....8.1385K]), yielding the precise spectral shape over four orders of
magnitude in rigidity [MalkMosk_2021]. To further test our model we turn to
other species.

In this letter, we argue that a single universal process is responsible for
the observed spectra of all CR species in the rigidity range below 100 TV. We
provide a simple formula, Eq. (4), that reproduces the spectra of all CR
species with only two shock parameters derived from a fit to the proton data.
We show how our formalism [MalkMosk_2021] operates using the proton, helium,
and carbon spectra and the B/C ratio.

Table 1: Model parameters and fit results for the proton spectrum. Parameter (St. err. %) |  $R_{0}$(GV) |  $R_{L}$(GV) | $q$ | $K=\left(\gamma+2\right)/\left(q-\gamma\right)$ |  $\chi_{\text{min}}^{2}$/dof | dof
---|---|---|---|---|---|---
Realistic Model (RM) | 5878 (3.5%) |  $2.24\times 10^{5}$ (28%) | 4.2 | 3.59 (4.9%) | 0.10 | 76-3
Loss-Free Model (LF) | 4795 (3.2%) | $\infty$ | 4.7 | 2.58 (2.9%) | 0.19 | 76-2
Table 2: Input parameters for CR species derived from their LIS [Boschini_2020]. Parameters | protons | helium | boron | carbon
---|---|---|---|---
$A_{s}$ (m-2 s-1 sr-1 GV-1) | $2.32\times 10^{4}$ | 3410 | 79 | 109
$\gamma_{s}$ | 2.85 | 2.76 | 3.1 | 2.76

## CR Reacceleration and Propagation Model.—

Only two combinations of parameters, constituting the bump amplitude, K, and
rigidity, $R_{0}$, suffice to predict the bump spectrum (formally requiring
six descriptors [MalkMosk_2021]). The reaccelerated CR density is elevated
over the background by the factor
$K=\left(\gamma+2\right)/\left(q-\gamma\right)$, where $\gamma$ is the
spectral index of the background CRs (normalized to $dR$) and $q>\gamma$ is
that of the shock, $q=\left(r+2\right)/\left(r-1\right)$. Here $r$ is the
shock compression ratio, which we determine from the proton fit. The enhanced
spectrum appears as a bump at some distance from the shock, as the
reaccelerated CR with lower rigidity do not reach the observer. There will
also be a cutoff to the CR enhancement, associated with the particle losses
from the flux tube and the limited shock size, which we discuss later.

The background CR spectrum is $f_{\infty}\propto R^{-\gamma},$ whereas the
shock would accelerate the freshly injected thermal particles to $R^{-q}$
spectrum. However, as the shock is likely to be weak, that is $q>\gamma\approx
2.85$ (for the protons), the freshly accelerated CRs are submerged into the
reaccelerated CR component. The CR density measured at a distance $z$ from the
shock is [BlandOst78, MalkMosk_2021]:

$\displaystyle f\left(R,z\right)$ $\displaystyle=$ $\displaystyle
f_{\infty}\left(R\right)\left[1+Ke^{-\Phi\left(R,z\right)}\right],$ (1)
$\displaystyle K$ $\displaystyle=$
$\displaystyle\left(\gamma+2\right)/\left(q-\gamma\right).$

The flux suppression exponent, $\Phi$, is a path integral of particles
propagating along the magnetic flux tube from the shock. One may distinguish
between two cases of the relative motion of the shock and the observer,
although the dependence of $\Phi\left(R\right)$ is similar in both cases. If
the shock normal makes an angle $\vartheta_{Bn}$ with respect to the local
magnetic field not too close to $\pi/2$, the CR distribution along the flux
tube is determined by the balance between the CR convection with the flow and
diffusion against it. In this case we have

$\Phi=u\int_{0}^{z}\frac{dz}{\kappa_{\parallel}\left(R,z\right)}+\mathcal{L}\left(R\right).$
(2)

Here $u=u_{sh}/\cos\vartheta_{Bn}$ with $u_{sh}$ being the shock velocity, and
$\kappa_{\parallel}$ is the parallel component of the diffusion tensor. The
term $\mathcal{L}\left(R\right)$ describes lateral losses from the flux tube.
The bump characteristic rigidity $R_{0}$ is introduced by the relation
$\Phi\approx$$\left(R_{0}/R\right)^{1/2}+\mathcal{L}$. Here the power $1/2$
corresponds to the Iroshnikov-Kraichnan turbulence derived as follows. The
parameter $R_{0}$ combines the distance to the shock, $z$, velocity $u$ and
the turbulence power, that defines $\kappa_{\parallel}$ in Eq. (2). According
to [MalkMosk_2021], $R_{0}\propto u^{2}z^{2}/l_{\perp}$ where the
characteristic scale of the CR pressure gradient across the flux tube is
$l_{\perp}^{-1}=P_{\text{CR}}^{-1}$$\partial P_{\text{CR}}$/$\partial r$. It
sets the amplitude of magnetic perturbations driven by this gradient.

In the special case $\vartheta_{Bn}\approx\pi/2$, we use the balance between
the divergence of convective and diffusive CR fluxes in the shock frame,
$\boldsymbol{v}\cdot\nabla
F\sim\nabla_{\parallel}\kappa_{\parallel}\nabla_{\parallel}F$, which is
consistent with Eqs. (1)–(2), and $F$ is the full-bump distribution related to
$f$ below. The CR diffusion proceeds primarily along the field but the
convective term is $\boldsymbol{v}\cdot\nabla
F\approx\boldsymbol{v}_{\perp}\cdot\nabla_{\perp}F$, which we may write
assuming the ISM wind being in $y$-direction and magnetic field in
$z$-direction. So we have v$\partial F/\partial
y\approx\kappa_{\parallel}\partial^{2}F/\partial z^{2}$, assuming that
$\kappa_{\parallel}$ is independent of $z$, for a simple illustration. The
result will then have a similar to the former case exponential dependence,
$\propto\exp\left[-vz^{2}/4y\kappa_{\parallel}\left(R\right)\right]$.
Averaging over $y$, or rather over radial coordinate inside the tube,
$r_{\perp}$, depends on the shock structure and CR propagation in the tube
that controls the losses. Thus, we consider the bump component
$f_{\text{b}}\left(R,z\right)\equiv f-f_{\infty}$ in Eq. (1) as a radially
averaged part of its full distribution, $F$. Namely,
$f_{\text{b}}=2a_{\perp}^{-2}\int_{0}^{a_{\perp}}F\left(R,z,r_{\perp}\right)r_{\perp}dr_{\perp}$,
where $a_{\perp}$ is the tube’s radius. The CR flux at $r_{\perp}=a_{\perp}$
is $-2\pi\kappa_{\perp}\partial F/\partial r_{\perp}|_{r_{\perp}=a_{\perp}}.$
Here $\kappa_{\perp}$ is a cross-field component of the CR diffusion tensor.
Calculation of the losses requires parameterization. Evaluating
$\kappa_{\perp}$ at $r_{\perp}=a_{\perp}$ is, perhaps, the most challenging
part. Nevertheless, we can express it through $\kappa_{\parallel}$ as
$\kappa_{\perp}\approx\kappa_{B}^{2}/\kappa_{\parallel}$, where
$\kappa_{B}=cr_{g}/3$ is the Bohm diffusion, with $r_{g}$ being the Larmor
radius. Using [MalkMosk_2021], one can represent the exponential decay of the
CR density along the tube as follows
$f\propto\exp\left[-\sqrt{R/R_{L}}\right]$, where $R_{L}$ is the rigidity cut-
off associated with the losses. It parameterizes the unknown flux through the
tube surface at $r_{\perp}=a_{\perp},$ and is roughly proportional to
$z^{-2}$. Therefore, the path integral in Eq. (2) at fixed $z$ can be
represented as

$\Phi=\left(R_{0}/R\right)^{1/2}+\left(R/R_{L}\right)^{1/2}.$ (3)

In this form, it makes the bump vanish also at high rigidities, $R\gg R_{L}$.
The loss-term is obtained by observing that $\kappa_{\perp}\propto R^{3/2},$
and estimating $\partial F/\partial r_{\perp}|_{a_{\perp}}\sim F/a_{\perp}$
for the CR flux across the magnetic tube boundary. The result in Eq. (3) then
follows from balancing this flux with the diffusive flux divergence in
$z$-direction. The parameter $R_{L}$ remains unknown in these simple
estimates, and we use it as the only adjustable parameter in the model, except
for $R\ll R_{L}$ where we can ignore the losses.

## Probing the Shock with Reaccelerated Protons.—

Following the formalism developed in [MalkMosk_2021], we represent the
spectrum of an arbitrary CR species in the following form, Eqs. (1), (3):

$f_{s}\left(R\right)=A_{s}R^{-\gamma_{s}}\left\\{1+\frac{\gamma_{s}+2}{q-\gamma_{s}}\exp\left[-\sqrt{\frac{R_{0}}{R}}-\sqrt{\frac{R}{R_{L}}}\right]\right\\},$
(4)

where $A_{s}$ and $\gamma_{s}$ are the normalization and spectral index of CR
species ($s=p$, He, B, C, …) that will be fixed using the data at $R\ll R_{0}$
(the so-called background CRs). These data are not related to the bump
phenomenon described by the remaining three parameters: $q$, $R_{0}$, and
$R_{L}$, where, $q=\left(r+2\right)/\left(r-1\right)>\gamma_{s}$.

In the first step we remotely sense the shock parameters using the precision
measurements of the local proton spectrum. The best available data in the
GV–TV domain are provided by AMS-02, CALET, and DAMPE. To mitigate the
systematic discrepancies between these instruments we use overlaps of their
rigidity ranges as follows. The most accurate measurements are provided by
AMS-02, which has several independent systems that allow for data cross
checks. To eliminate the solar modulation effects below the first break
rigidity, we use the proton LIS spectrum [Boschini_2020], which was derived
using the AMS-02 data. The CALET and DAMPE data are then modified to normalize
their first break rigidity $R_{{\rm br}}^{\prime}$ and their flux at $R_{{\rm
br}}^{\prime}$ (minimum in the spectrum in Fig. 1) to the values derived from
the proton LIS. This was achieved by applying the following transformations:
$R=\tilde{R}/1.12$, Flux$=\widetilde{{\rm Flux}}\times$1.09 to CALET, and
$R=\tilde{R}/0.95$, Flux$=\widetilde{{\rm Flux}}\times$0.98 to DAMPE data.
Here $\tilde{R}$ and $\widetilde{{\rm Flux}}$ are the published rigidity and
flux values [2019PhRvL.122r1102A, eaax3793]. These adjustments are well within
the instrumental error bars. The model parameters derived from the best fit to
the modified data of all three instruments are summarized in Table 1.

The LIS spectra below the break were used to obtain the model input parameters
listed in Table 2. Fig. 1 shows the proton data along with two models, where
we have considered the cases of $R_{L}=\infty$ (LF model) and $R_{L}<\infty$
(RM model), more realistic for $R\gtrsim R_{{\rm br}}^{\prime\prime}$. The RM
model provides a better match of the data above 10 TV indicating that the
shock size is limited. We emphasize that all model parameters are now fixed.

Figure 1: Proton rigidity spectrum from AMS-02, CALET, and DAMPE [Aguilar2021,
2019PhRvL.122r1102A, eaax3793]. CALET and DAMPE data were modified to adjust
for systematic differences between the three instruments, see text for
details. The curves show the model calculations, Eq. (4), for the cases of
$R_{L}=\infty$ (LF model) and a more realistic for high rigidities case of
finite losses $R_{L}<\infty$ (RM model). The model and input parameters are
summarized in Tables 1, 2.

## Calculating Spectra of Other Elements.—

According to Eq. (4), the species are differentiated solely by their spectral
index $\gamma_{s}$ below the break $R_{{\rm br}}^{\prime}$. We illustrate this
in Fig. 2, where we use the spectral indices corresponding to the helium and
boron LIS, and an arbitrary index corresponding to a steeper spectrum,
$\gamma_{{\rm arb}}=3.7$. We see that CRs with steeper background spectra
undergo stronger reacceleration.

Fig. 3 shows the primary helium and carbon spectra. The DAMPE and CREAM helium
spectra were normalized to the helium LIS by applying factors 1.08 and 1.05,
correspondingly. Following [2020PhRvL.125y1102A] the CALET spectrum of carbon
was multiplied by 1.27 to agree with AMS-02 data. Given the same spectral
indices below the first break $R_{{\rm br}}^{\prime}$, the spectral shapes are
identical.

We are also comparing our model calculations with the measurements of the B/C
ratio that is widely used to evaluate the CR residence time in the Galaxy. It
is measured in a wider rigidity range and with higher precision than the
individual spectra of boron and carbon. Fig. 4 shows the calculated B/C ratio
compared with available data. The inset shows the ratio multiplied by a factor
$R^{0.3}$ to emphasize the fine details in the bump range. Our model shows
$\approx$10% deviation from predictions for the Kolmogorov turbulence
[1941DoSSR..30..301K] ($\propto R^{-0.33}$).

Figure 2: The bump strengths in the CR spectra with different spectral indices
corresponding to the helium and boron LIS (Table 2), and $\gamma_{{\rm
arb}}=3.7$ (for illustration), calculated using Eq. (4) and compared to the
straight background power-laws. The ratios of the enhanced fluxes at $R=30$ TV
to the background values for the three different values of $\gamma$ are 2.46,
3.06, and 6.05, respectively.  Figure 3: The calculated helium and carbon
spectra are compared with available helium data from AMS-02, DAMPE, CREAM,
ATIC, and NUCLEON and carbon data from AMS-02, CALET (multiplied by 1.27),
ATIC, CREAM, and NUCLEON [Aguilar2021, 2020PhRvL.125y1102A,
2020PhLB..81135851K, 2017ApJ...839....5Y, 2009BRASP..73..564P,
2020JETPL.111..363K, 2009ApJ...707..593A, DAMPE_He]. The RM model calculations
are using Eq. (4) with parameters provided in Tables 1, 2.  Figure 4: The
calculated B/C ratio is compared with AMS-02, CREAM, ATIC, and NUCLEON data
[2019AdSpR..64.2559G, Aguilar2021, 2008APh....30..133A, 2008ICRC....2....3P].
The RM model calculations are using Eq. (4) with parameters provided in Tables
1, 2.

## Background spectra.—

Here we comment on the origin of the background spectra of different CR
species below $R_{{\rm br}}^{\prime}$. A physical explanation for the observed
difference between $\gamma_{p}$ and $\gamma_{s}$ of primary elements with
effective mass to ion-charge ratios $\left\langle A/Q\right\rangle$ $\approx
2$, or higher has already been suggested [m98, MDSPamela12, Hanusch2019ApJ].
They develop harder than proton spectra (the index difference
$\gamma_{p}-\gamma_{s}\approx 0.1$) because SNR shocks extract them from
thermal plasma most efficiently when these shocks are at their strongest. The
time-integrated spectra of heavier elements become then harder than that of
the protons.

This bias stems from the shock related turbulence whose amplitude increases
with the Mach number. Resonantly driven by protons, it curbs them tighter
downstream than heavier nuclei, so the latter may cross the shock to gain
energy. As the above differentiation mechanism depends exclusively on
$\left\langle A/Q\right\rangle$, its salient consequence is that all elements
with the same $\left\langle A/Q\right\rangle$ must have the same
ultrarelativistic rigidity spectra [Malkov2018a]. No exceptions to this rule
have been reported so far [Aguilar2021].

The steep background spectrum of the secondaries is related to a decreasing
residence time of CR species in the Galaxy as rigidity increases. The
secondaries are already produced in fragmentations of primary CRs with a steep
spectrum that is further steepened due to the rigidity-dependent residence
time.

## Discussion.—

The observed bump in the CR spectra is not unique. Similar bumps may be
generated in other locations where shocks are present. Whether their integral
contribution competes favorably with the CR reacceleration by the ISM
turbulence is too early to say. Theoretical estimates of the turbulent
reacceleration [1994ApJ...431..705S, 2006ApJ...642..902P, DruryStrong2017]
suggest that it supplies up to 50% of the power required to maintain the CRs
against losses, but it is hard to prove on observational grounds, as these CRs
are inseparable from the background.

By contrast, a contribution of CR-reaccelerating shock is identifiable and, as
we argued, is currently observed. The next step is then to create an inventory
of such shocks and integrate their contributions. Our results show that such
shocks are able to temporarily modify the local CR spectrum. This is not
entirely new, e.g., see a dedicated study of time-dependent CR transport
[2019ApJ...887..250P]. However, a truly unique result is that we do not need
such a major event as a supernova explosion nearby to modify the local CR
spectrum.

Interestingly, the B/C ratio is also affected by such modifications. One can
thus imagine that the observed B/C ratio cannot be used for estimates of the
CR confinement time without appropriate corrections for the reacceleration
effects due to the passing stars. We can further speculate that even the well-
known “knee” at $\sim$3 PV could be the reacceleration enhancement due to a
larger by 2-3 orders of magnitude shock.

###### Acknowledgements.

We thank Samuel Ting, Andrei Kounine, and Vitaly Choutko for discussions of
the AMS-02 data, John Krizmanic, Kazuyoshi Kobayashi, and Shoji Torii for
discussions of the CALET data, Jin Chang and Qiang Yuan for providing the
DAMPE data prior to their publication and discussions, Dmitry Podorozhny,
Alexander Panov, and Andrey Turundaevsky for providing the updated data of
NUCLEON and ATIC, and Eun-Suk Seo for providing the CREAM data. Mikhail A.
Malkov acknowledges support from NASA ATP-program within grant 80NSSC17K0255
and from the National Science Foundation under grant No. NSF PHY-1748958. Igor
V. Moskalenko acknowledges support from NASA Grant No. NNX17AB48G.


#  The Eight Epochs of Math as regards
past and future Matrix Computations

Frank Uhlig
Department of Mathematics and Statistics,
Auburn University, AL 36849-5310<EMAIL_ADDRESS>
Abstract :  This talk gives a personal assessment of Epoch making advances in
Matrix Computations from antiquity and with an eye towards tomorrow.
We trace the development of number systems and elementary algebra, and the
uses of Gaussian Elimination methods from around 2000 BC on to current real-
time Neural Network computations to solve time-varying linear equations.
We include relevant advances from China from the 3rd century AD on, and from
India and Persia in the 9th century and discuss the conceptual genesis of
vectors and matrices in central Europe and Japan in the 14th through 17th
centuries AD.
Followed by the 150 year cul-de-sac of polynomial root finder research for
matrix eigenvalues, as well as the superbly useful matrix iterative methods
and Francis’ eigenvalue Algorithm from last century.
Then we explain the recent use of initial value problem solvers to master
time-varying linear and nonlinear matrix equations via Neural Networks.
We end with a short outlook upon new hardware schemes with multilevel
processors that go beyond the 0-1 base 2 framework which all of our past and
current electronic computers have been using.

## Introduction

In this paper we try to outline the Epoch making achievements and
transformations that have occurred over time for computations and more
specifically for matrix computations. We will trace how our linear algebraic
concepts and matrix computations have progressed from the beginning of
recorded time until today and how they will likely progress into the future.
We take this limited tack simply because in modern times matrices have become
the elemental and universal tools for most any computation.
This evolution of our matrix methods will be described in broad strokes. My
main emphasis is to trace the mathematical genesis of Matrices and their uses
and to learn how the modern matrix concept has evolved in the past and how it
is evolving. I am not interested in Matrix Theory by itself, but rather in
Matrix Computations, i.e., how matrix concepts and algorithms have been
developed from approximately 3000 BC to today, and even tomorrow.
This paper describes eight noticeably separate Epochs that are distinguished
from each other by the introduction of evolutionary new concepts and
subsequent radically new computational methods. Following the historical trail
through six historically established Epochs, we will then look into the
present and the near future.
What drives us to conceptualize and compute differently now, what is leading
us into the 7th and possibly 8th Epoch? When and how will we likely compute in
the future.
I am not a math historian, I have never taught a class in math history.
Instead throughout my academic career I have worked with matrices: in matrix
theory, in applications and in numerical analysis and I like to construct
efficient new algorithms that solve matrix equations. The idea for this paper
is in part due to my listening by chance to a very short English broadcast
from Egyptian Radio on short wave some 40 years ago in the 1970s. It described
an Egyptian papyrus from around 2000 BC that dealt with solving linear
equations by row reduction and zeroing out coefficients in systems of linear
equations, i.e., by what we now call “ Gaussian Elimination”. When I heard
this as a young Ph. D., I was fascinated and wrote the station for more
information. They never answered and when I was in Cairo many years later, the
Egyptian Museum personnel could not help me either with locating the source.

Thus I became aware that Carl Friedrich Gauß did not invent what we now call
by his name; but who did?
For many decades, this snippet of math history just lingered in my mind until
a year ago when I was sent a book on Neural Network (NN) methods for solving
time-varying linear and non-linear equations and was asked to review it. The
NN methods were - to me and my understandings of numerics then - so other-
worldly and brilliant that I began to think of the incredible leaps and
’bounces’ that math computations have gone through over the eons, from era to
era. I eventually began to detect 7 or 8 computational Sea Changes, what I
call “ Epochs”, in our ability to compute with matrices and that is my topic.

##  1 A Short History of Matrix Computations

Nobody knows how numbers and number systems came about, just like nobody knows
‘who invented the wheel’. I will start with a few historical facts about
number systems and how they developed and were used across the globe in
antiquity.

### 1.0 Early Number Systems

Humankind’s first developments of number systems were very diverse and
geographically widely dispersed, yet rather slow. The first circle cypher for
zero occured in Babylonia around 2500 BC, or 4500 years ago. A continent or
two removed, the Mayans used the same concept and circle zero symbol from
around 40 BC. In India it was recognized during the 7th century. But zero only
became recognized as a ‘number to compute with’ like all the others in the 9th
century in central India. Our decimal system builds on the ten numbers 0, 1,
2, 3, 4, 5, 6, 7, 8, and 9. The decimal positional system came from China via
the Indus valley and it started to be used in Persia in the 9th century. It
was combined with or derived from a Hindu number system of the same time
period.
In fact westerners call the current decimal number symbols wrongfully
’arabic’, but most westerners (and I) cannot read the license plates on cars
in Egypt since the arabic world does not use our Persian/Hindu numbers in
writing but its own script using arabic letters to designate numbers. Should
we call our ’western’ numbers ’farsi’ or ‘hindu’ instead?

Various bases have been used for numbering. There have been base 2, base 8,
base 10, base 12, base 16, base 60, base 200 number systems and possibly more
at some time somewhere throughout human history. Counting and simple
computations started with notched sticks for record keeping and with the
invention of sand or wax tablets and then the abacus. These simple tools were
developed a little bit differently and independently in many parts of the
globe.

###  1.1 Antiquity (1st Epoch)

Around 2200 to 1600 BC, Sumerian, Babylonian and Egyptian land survey
computations became mathematized in order to mark and allot land after the
yearly Euphrates, Tigris and Nile floods. That naturally lead to linear
equations in 2, 3 or 4 variables and subsequent methods to solve them that
amounted to what we now call row reduction or Gaussian elimination.
Mathematics Computations did not advance much during the Greek times as Greek
mathematicians were mainly interested in mathematical theory and in
establishing the concept of a formal proof, as well as elementary number
theory of which the Euclidean algorithm is still used today.
Neither did the complicated Roman numerals lend themselves to easy
computations and no further computational advances happened there.

###  1.2 Early Mathematical Arts in China, India and the Near East (2nd Epoch)

(Based in part on a lecture at Hong Kong University in 2017, given by Zhou
Xiangyu, Chinese Academy of Sciences, for Chinese sources. And on Indian and
Arabic sources from elsewhere)
In prehistoric and historic times (1600 BC - 1400 AD), knot and rod calculus
were prevalent in China. They were based on a decimal positional system, the
so-called ”rod numerals”. These comprised the most advanced number system of
the time and it was used for several millennia before being adopted and
expanded in Persia and India in the 9th century AD and later on adopted in
central Europe.
The ”Mathematical Classic of Sunzi” by Suanjing (from the 3rd - 5th century)
gives a detailed description of the arithmetic rules for counting rods. In the
Indus valley clay tablets covered with sand were used for mathematical
computations several millennia ago. Bhaskara (600 - 680 AD) in India was first
to write numbers in the Hindu positional decimal system which used the circle
for zero. In 629 AD he approximated the sine function by rational expressions
while commenting on Aryabhatta’s (476 - 550 AD) book ”Aryabhatiyabhisya” from
499 AD. An Indian contemporary of Bhaskara, Brahmagupta (598 - 665 AD) was the
first to establish the rules that govern computing with zero. Brahmagupta
texts were written in Sanskrit verse that used the Sanskrit word for ’eyes’ to
denote 2 and ’senses’ for the number 5 etc. This was common in Indian
mathematics and science writings at the time. The earliest record of
multiplication and division algorithms using the Hindu numerals 1 through 9
and 0 was in writings by Al Khwarizmi 780 - 850 AD, a Persian mathematician
employed in Bagdad. His ”Book of Manipulation and Restoration” established the
golden rule of Algebra that an equation remains true if one subtracts the same
quantity from both sides. He also wrote down multiplication and division rules
that are identical to those of Suanjing from the 3rd to 5th century and China.
To Suanjing we also owe the Chinese Remainder Theorem. Finally the advanced
Hindu-Arabic decimal number system was introduced into the West by Leonardo
Fibonacci (1175 - 1250 AD) of Pisa in his ”Liber Abaci” or ’The Book of
Calculations’ (1202),
Applied and numerical computations were driving much of Chinese mathematics.
Wang Xiaotang (580 - 640 AD), for example, tried to find the roots of cubic
polynomials that appeared in civil engineering and water conservation
problems. In the Mathematical Treatise in Nine Sections of 1247, Qin Jiushao
(1202 - 1261) developed the ’Qin Jiushao method’ which is now commonly called
the ’Horner-Ruffini scheme’ for computing with and finding roots of
polynomials iteratively. William George Horner (1819) [17] and Paolo Ruffini
(1804-1807-1813) reinvented the Qin Jiushao method unknowingly six hundred
years later.
Chinese rod calculus was the method of choice for computing in China until the
abacus took over during the Ming Dynasty 1388 - 1644. Cheng Dawei 1355 - 1606
is the author of the first ’Numerical Analysis’ book titled ”The General
Source of Computational Methods” and published in 1592. It describes methods
to add, subtract, multiply, and divide on an abacus. The abacus itself was
invented in various incarnations at various times and in several locations of
the globe. It essentially combines several decimal rods on one board with
beads on strings.
Chinese mathematicians from the 3rd century BC onwards to the 10th century AD
brought us the ”Nine Chapters on the Mathematical Arts” that uses the numbers
1 through 9. This book was later disseminated further to the west, to India
and Persia as described above. In its 7th chapter, determinants first appeared
conceptually, while chapter 8 abstracts the concept of linear equations to
represent them by matrix-like coefficient tableaux. These ‘matrix equations’
were solved in China, again by ‘Gaussian elimination’, 1500 years before Gauß’
birth and 1800 years after the middle-eastern seasonal flood prone countries
had first used the Gaussian algorithm around 1800 BC. Gauß himself described
the method as the “common method of elimination” in his papers and
mathematicians then attached his name to it as an honor.

###  1.3 The Genesis of Vectors and Matrices (3rd Epoch)

To advance matrix computations further there was a need to conceptualize
coordinates and vectors in space.
In the 14th century AD, Nicole Oresme developed a system of orthogonal
coordinates for describing Euclidean space. This idea was taken up by René
Descartes in the 17th century and is familiar to all of us now under the
concept of Cartesian Coordinates. Thereby the world became ready for matrices
and matrix computations in their own right.
In 1683 Gottfried Leibnitz in Germany and Seki Kowa in Japan both and
unbeknownst to each other re-invented the concept of a matrix as a rectangular
array of coefficients for studying linear equations. Leibnitz also used and
suggested row elimination to simplify and find their solutions. These efforts
enabled Gauß to repeat what the Egyptians had done four millennia earlier: He
was asked to survey the lands of his ruler, the Archduke George Augustus of
Hanover and measure the size of this kingdom inside Germany in the early
1800s. Beginning in 1821 and 1823 Gauß , as Professor of Geodesy (and not of
Mathematics) in Göttingen, would measure the angles and distances between many
of the highest points there, such as the Brocken, the Inselsberg, 104 km
apart, and the hills around Göttingen and later expanded the surveys all the
way to the North Sea. He did this multiple times, preferably when the weather
was clear. Thereby he set up systems of linear equations with generally more
equations than unknown due to repeated measurements on different days.
To solve these overdetermined and naturally ‘unsolvable’ system $Ax=b$, Gauß
devised the normal equation $A^{T}Ax=A^{T}b$ and solved it approximately. But
the normal equations method eventually turned out to be numerically unsound.
It took over a century to find out why, the reason being that condition
numbers multiply, see Olga Taussky (1950) [27].

###  1.4 Eigenvalues and the Characteristic Polynomial (4th Epoch)

As differential operators and matrices were beginning to be investigated and
dealt with by the early 1800s, their connections and similarities were slowly
recognized in the mathematics world.
The replication of certain functions $f\neq 0$ by a given differential
operator $\cal A$ was noticed first and became the subject of studies. What
were the functions $f$ for which ${\cal A}f=\alpha f$ for some scalar
$\alpha$? How could they be found from $\cal A$, what about $\alpha$?
In 1829 Augustin Cauchy [6, p. 175] began to view the erstwhile ‘eigenvalue
equation’ ${\cal A}f=\alpha f$ as a ‘null space equation’, namely ${\cal
A}f-\alpha f=0$ or $({\cal A}-\alpha\ id)f=0$ for the identity operator with
$id\ f=f$ for all $f$.Complete knowledge of the eigenvalues $\alpha$ and
eigenfunctions $f$ of a differential operator $\cal A$ allowed for a simple
sum representations of the general solution of the linear differential
equation described by $\cal A$. Thus Cauchy’s ‘null space equation’ became
essential for determining the behavior of systems governed by linear
differential equations.
Cauchy’s knowledge of and interest in determinants (think Cauchy-Binet
Theorem) then led him to define the ‘characteristic polynomial’ of a square
matrix $A$ in 1839 [7, p. 827] as $f_{A}(\alpha)=\det(A-\alpha I)$ and thereby
he initiated renewed studies in polynomial root finding algorithms in the hope
of obtaining analogous diagonalization results for linear matrix times vector
products. And the search for polynomial rootfinders was on. By modern day
hindsight, reducing the eigenvalue problem from an $n^{2}$ data problem of the
entries of a matrix $A_{n,n}$ to one of the $n+1$ coefficients of its
characteristic polynomial is data compression and therefore it was doomed to
fail. But that remained unrealized by the mathematics community for more than
100 years.
James Sylvester finally gave the tableau concept of matrices its name ’matrix’
in 1848 or 1850. And after roughly 2 decades _1829 – $>$ 1839 –$>$ 1848/50_,
the first century of Matrix Theory or theoretical Linear Algebra had begun.

Back to matrix computations:
Cauchy’s idea led mathematicians to try and compute the characteristic
polynomials of matrices and find their roots in order to understand the eigen-
behavior of matrices. We still teach many concepts and lessons today that are
based on the ’characteristic polynomial’ $f_{A}(x)$ of a matrix $A$. Why, we
should ask ourselves. Because unfortunately studying ’characteristic
polynomials’ in place of matrices has turned out be a costly dead end for
computational and applied mathematics: In the century and a half that followed
Cauchy’s work, more than 4,000 papers on computing the roots of polynomials
were published, together with 200 to 300 books on the subject, bringing us
many algorithms, all of which failed more often than not. Many illustrious
careers and schools of mathematics were founded based on this unfortunate and
ever elusive goal.
During the same period, 2-D hand crank computing machines were invented and
built to effect long number multiplications and divisions. First by Charles
Babbage, then as commercial geared adding machines that were still being used
in office work well into the 1960s. These worked as two-dimensional abaci of
sorts. But eventually digital (at first punch-card fed) computers became the
tools of our computational trade in science, in engineering, in business, in
GPS, in Google, in social media, large data, automation, etc, etc.

But how could we or would we find matrix eigenvalues accurately? A turn-
around, a new method, a new computational Epoch was needed. From where, by
whom, and how …?

###  1.5 Iterative Matrix Algorithms (5th Epoch)

To move us forward, it appears that matrix methods themselves might have to be
developed that would solve the matrix-intrinsic eigenvalue problem by
themselves. But before that was possible there were further unfortunate
‘detours’.
Carl Friedrich Gauß - in his doctoral thesis of 1799 [13] \- had disproved all
earlier attempts to establish the Fundamental Theorem of Algebra, i.e., that
all polynomials over the reals numbers can be factored into as many real or
complex conjugate factors as their degree says. His thesis then included the
first complete and correct proof of the Fundamental Theorem of Algebra.
In 1824 Niels Abel [2] showed that the roots of some 5th degree polynomials
can not be found by using radical expressions of their coefficients; Gauss
never opened or read the submitted paper and thus in fact rejected it
knowingly on the ‘grounds’ that God would not have complicated the World thus
… for us. Abel published his result privately, a broken man. Évariste Galois
[10, 11] extended Abel’s result in 1830 by giving group theoretic conditions
for polynomials to be solvable by radicals; the extended paper (introducing
Galois Theory) was originally rejected and appeared only posthumously in 1846
[12].
Cautioned by these ‘rejected’ inconvenient results, the polynomial approach to
matrix eigenvalue computations could have been shunned by clearer heads early
on; but the ‘dead end’ determinants and characteristic polynomial roots road
was taken instead for more than a century. Note that Cauchy’s matrix result
and most other fundamental matrix results from the 19th century were
formulated in terms of determinants and only in the mid 20th century did the
term ’matrix’ appear in matrix theoretical article and book titles.
A matrix based approach to the eigenvalue problem nowadays starts from the
simple fact that for any $n$ by $n$ matrix $A$ and any $n$ vector $b$, the
sequence of vector iterates $b,Ab,A^{2}b,...,A^{k}b,...,A^{n}b$ contains $n+1$
vectors in $n$-space which makes these vectors linearly dependent. Their
linear dependency then leads to an $n$th degree polynomial $p_{b}(A)$ that
sends $b$ to zero. The vanishing polynomial for any $b$ turns out to always be
a factor of the characteristic polynomial of $A$ and it can be found by
Gaussian elimination rather than using determinants.
The same idea shows that vector iteration converges for every starting vector
$b\neq 0$ and any given square matrix $A$ and this has lead engineers in the
early 20th century to construct iterative matrix algorithms that could solve
linear equations and the matrix eigenvalue problem. Iterative matrix
algorithms actually do go back further, such as to the Jacobi (1839) [18],
Gauß -Seidel (1874) , and various SOR methods that are designed to solve
linear systems iteratively. These use matrix splittings of $A$ rather than
vector iteration.
Alexei Krylov (1931) [20] introduced and studied the vector iteration
subspaces span$\\{b,Ab,...,A^{k}b\\}$ in their own right. Following his ideas,
large sparse matrix systems are nowadays treated iteratively in so called
Krylov-based methods, both to solve linear equations and to find matrix
eigenvalues. Standard widely used Krylov type iterative matrix algorithms
carry the names of Steepest Descent, Conjugate Gradient, Arnoldi (1951) [3],
Lanczos (1950) [21]. Others are called GMRES, BICGSTABLE, QMR, ADI etc. Most
Krylov type methods are matrix and problem specific and they are mostly used
for huge sparse and structured matrices where direct or semi-direct methods
cannot be employed due to their high computational and storage costs. Krylov
methods generally rely on preconditioners $M$ for linear systems $Ax=b$ that
shift the spectrum of $M^{-1}A$ for faster convergence and they thrive on
incomplete matrix splittings etc. Typically they give only partial results.
Who needs to know all the million eigenvalues of a million by million matrix
model. Krylov methods can be tuned to give information where it is needed for
the underlying system.

###  1.6 Francis Algorithm and Matrix Eigenvalues (6th Epoch)

The Second World War and post Second World War periods were filled with
innovations:
The atomic era had begun, as well as rocket science; commercial air flight
became popular; digital computers were being developed, first as valve
machines and later transistorized. Supersonic speeds were realized, Computer
Science was developed, etc. But there were many crashes and disasters with the
new technologies: Commercial aircraft (Super-Constellation, Convair, …) and
military ones (Starfighter …) would crash weekly around the globe; and newly
built suspension bridges would collapse in strong winds.
The crux of the matter was that while matrix models of the underlying
mechanical systems could readily be made using the laws of physics and
mechanics, no one could reliably compute their eigenvalues. Engineers could
not test their designs for eigen-modes in the right half plane! And Krylov
methods were unfortunately not sufficient for testing for eigenvalues in a
half plane.
If a matrix model of a mechanical or electrical or … structure has right-half
plane eigenvalues $\lambda$ then - upon proper excitation - there would be a
eigen component of the ever increasing form $e^{\lambda t}\to\infty$ as
$t\to\infty$ that resonates and self-amplifies inside the structure itself.
This then leads to ever increasing destructive vibrations and ultimate
failure. The aircraft ’flutter problem’ was discovered during World War 2. In
England during WW2, Gershgorin circles that contain all of a system’s
eigenvalues were drawn out in the complex plane by rather primitive valve
computers and checked to ascertain system stability.
The general matrix eigenvalue problem was finally solved independently and
similarly by John Francis in London and by Vera Kublanovskaya in Russia nearly
simultaneously around 1960. Francis’ (or the QR) algorithm [8, 9] is based on
Alston Householder’s idea to try and solve matrix problems by matrix
factorizations. Francis’ method is an orthogonal subspace projection method
and it works differently than the Krylov based methods which solve a given
matrix eigenvalue problem by projecting onto a Krylov subspace that is derived
from and suitable for $A$.
A ’divide and conquer’ matrix factorization strategy was first employed by
Heinz Rutishauser (1955, 1958) [25, 26] in his LR matrix eigenvalue algorithm:
If one can factor $A=LR$ into the product of a lower and an upper triangular
matrix $L$ and $R$ as $A=LR$ and if $L$ is invertible, then for the reverse
order product $A_{1}=RL$ we have $A_{1}=L^{-1}AL$ since $R=L^{-1}A$. If
$A_{1}$ again allows an LR factorization $A_{1}=L_{1}R_{1}$ with $L_{1}$
nonsingular, then by reverse order multiplying we obtain

$A_{2}=L_{1}^{-1}A_{1}L_{1}=L_{1}^{-1}L^{-1}ALL_{1}$

and so for the sequence of likewise constructed matrices $A_{i}$ for $i=3,...$
if the respective LR factorizations are possible at each stage $i$. In this
case the iterates $A_{i}$ clearly remain similar to the original matrix $A$
and thus the iterates all have the same eigenvalues as $A$. Note, however,
that if for example the (1,1) entry $A(1,1)$ is zero in $A$, then there exists
no un-pivoted LR factorization for A and the method breaks down since pivoting
is a one sided matrix process and not a similarity. Therefore Rutishauser’s
method is only applicable to a limited set of matrices $A$ for which every LR
iterate $A_{i}$ allows an un-pivoted LR factorization. Rutishauser had noted
that if LR factorizations are possible for all iterates $A_{i}$, then the
$A_{i}$ become nearly upper triangular for large $i$ with $A$’s eigenvalues on
the diagonal. (Very loosely said.)
John Francis was very interested in the flutter problem at the time when, by
chance, someone dropped Rutishauser’s 1958 LR paper [26] on his desk at the
CRDC in London.  _(In my interview with John Francis in 2009[16], he did not
know who that might have been.)_ Francis was aware through contacts with Jim
Wilkinson of the backward stability of algorithms that involve orthogonal
matrices Q. So rather than using Gaussian elimination matrices $L$, Francis
experimented with orthogonal $A=QR$ factorizations. At roughly the same time
Vera Kublanovskaya worked on an LQ factorization of $A$ as $A=LQ$ and
subsequent reverse order multiplications [19] in Leningrad, Russia, that would
also compute the eigenvalues of $A$.
Rutishauser had observed convergence speed-up for his LR method when replacing
$A$ by $A-\alpha I$, i.e., shifting. Hence Francis experimented with shifts
for QR and then established the ’Implicit Q Theorem’ [9] in order to
circumvent computing eigenvalues of real matrices over the complex numbers.
Implicit shifts also avoid rounding errors that would be introduced by
explicit shifts. Francis’ second paper (1962) [9] also contains a fully
computed flutter matrix problem of size 24 by 24. The eigenvalues of such
’large’ problems had never before been computed successfully.
Francis’ Implicit Q Theorem then allowed Gene Golub and Velvel Kahan (1965)
[15] to compute singular values of large matrices for the first time and this
application later spawned the original Google search engine and brought us -
in a way - into the internet age.
In 2002 the Multishift QR Algorithm was developed by Karen Braman, Ralph
Byers, and Roy Mathias [4, 5]. It relies on subspace iteration and extends
Francis’ QR and combes it with Krylov like methods. This extension allows us
today to compute the complete eigen and singular value structure of dense
matrices of sizes up to 10,000 by 10,000 economically on laptops.

What is being missed today computationally? What Epoch(s) might come next? Why
and how?

###  1.7 Two New Epochs Ahead (7th and 8th Epochs) – yet to come

Two new Epoch generatinging impulses have become visible on the matrix
computational horizon of today:
(A)  One expands our computational abilities from static problem solving
algorithms to real time methods for time-varying problems, and
(B)  the other involves computer hardware advances.

1.7.1 Time-varying Problems and Real-time Solvers  _( Epoch 7)_ :
Our current best numerical codes can solve static problem very well; that is
what they are designed for.
As we begin to rely more and more on time-dependent sensoring and on robotic
manufacture, we need to learn how to solve our erstwhile static equations, but
now in real time and with time-varying coefficients, and preferably accurately
as well. It seems quite alluring to try and solve a time-varying problem by
using the static time-dependent inputs at each instance statically. But such a
naive solution cannot suffice since at the next time step, whose ‘solution’ we
have just computed ‘statically’, the problem parameters have already changed
and thus our ‘static’ solution solves a completely different problem, which –
unfortunately – has little value. _If any at all_.

1.7.2 Computer Hardware  _( Epoch 8 )_ :
Since the earliest electronic computing devices of the 1940s, all our
computers have worked as giant and embellished Turing machines with logic
gates, switches and memory that rely only two numerical states: 0 and 1 or on
or off. Hence all our computer data is stored and manipulated as sequences of
0 and 1.
Lately our computing ability has come up against the limits of storing and
working with data and processors that can only deal with zeros and ones. Our
processing speeds have not advanced significantly over the last couple of
years; we are still stuck with 3 - 4 GHz processors. To alleviate this
bottleneck, chip makers have created multi-processor chips and software firms
have introduced better and quicker software and operating systems, but the
basic processor speed has not budged much.
At this time computer scientist and manufacturers are trying to overcome this
0 - 1 bottle-neck by replacing our 0 - 1 processors, chips, memories and
transistors by improved transistors and chips that can store and process
multi-states, such as 0-1-2-3-4 or 0-1-2-3-4-5-6-7-8 or even higher numbered
data representations. This could lead us to another computing sea change
bringing us into a new computational Epoch via hardware. And further out on
the horizon lies the possibility of having infinitely many quantum states
based computers.

##  2 On Neural Network Methods _(Epoch 7_ already under way)

The last century brought us valuable tools to solve most static problems that
involve matrices.
Our current numerical matrix tools can solve static linear equations, matrix
equations such as Sylvester or Lyapunov equations, eigenvalue problems and
generalizations of these, both of the dense or structured and of the solvable
or unsolvable kind; likewise we can solve static optimization problems of all
sizes and for nearly all structured matrices and for most static applications.
But what do do with such problems when their coefficients are time-varying or
time dependent?
In numerical computations, there has always been a see-saw between models that
resulted in derivative inspired differential equations and in linear algebra
based matrix equations. Their respective computational advantages differed
from problem to problem. Neural networks (NN) are an amalgam of matrix methods
and differential methods using a mixture of both. NN methods are designed to
solve time-varying dynamical systems. Numerical methods for time-varying
dynamical systems first came about in the 1950s and subsequently have gained
strength in the 1980s and 1990s and beyond, see the introduction in Getz and
Marsden [14] for example. To solve a time varying equation $f(X(t))=G(t)$,
these studies start from the error equation $E(t)=f(X(t))-G(t)$ and stipulate
exponential decay of the error function $E(t)$ by trying to solve

$\dot{E}(t)=-\lambda E(t)$ (1)

for a positive decay constant $\lambda$. There are essentially three ways to
go about solving the error function differential equation (1): homotopy
methods, gradient methods and ZNN Neural Network methods.

###  2.1 A Neural Network approach to solve time-varying linear equations
$\mathbf{A(t)x(t)=b(t)}$

Here $A(t)$ is a nonsingular time-varying $n$ by $n$ matrix and
$b(t)\in\mathbb{R}^{n}$ is a time-varying vector, respectively. Clearly the
unknown solution $x(t)$ of the associated linear equation $A(t)x(t)=b(t)$ will
be time-dependent as well.
The first paper on Zhang Neural Networks (ZNN) was written by Yunong Zhang et
al. in 2002, see [28]. Today there are over 250 papers, mostly in engineering
journals that deal with time-varying applications of the ZNN method, either in
hardware chip design for specialized computational tasks as part of a plant or
machine or for time-varying simulation problems in computer algorithms and
codes. Unfortunately, the ZNN method and the ideas behind ZNN are hardly known
today among numerical analysis experts. The method itself starts with using
Suanjing’s and Al Khwarizmi’s rule for reducing equations that first appeared
1 1/2 millennia ago. This simple rule was also employed by Cauchy in 1829 to
transform the static matrix eigenvalue problem from $Ax=\lambda x$ to
$Ax-\lambda x=0$ and finally to det$(A-\lambda I)=0$. For the time-varying
linear equations problem Zhang’s Neural Network method starts with

$A(t)x(t)-b(t)=0\ .$

ZNN then looks at the error function

$E(t)=A(t)x(t)-b(t):\mathbb{R}\to\mathbb{R}^{n}$

Note that standard static methods would look at an error norm $\|E\|$ for the
error function. Neural Networks do not, they study the error function $E(t)$
instead. And they start with an implicit ‘ideal wish’: _What could or should
we wish for $E$_? Time-varying computations would be near ideal if their error
functions were decaying exponentially fast. This is impossible to achieve (or
even ask for) with our best static equations and problem solvers of the 21st
century. For static numerical matrix methods backward stability is considered
best.
In all Zhang NN methods we stipulate that the error function $E(t)$ decreases
exponentially fast over time to the zero function. This means that

${\dot{E}(t)=-\lambda E(t)}\ \text{ for some chosen number }\lambda\gg 0,\
\text{ the decay constant}.$

Note that for the time-varying linear equations problem we have

$\dot{E}(t)=\dot{A}(t)x(t)+A(t)\dot{x}(t)-\dot{b}(t)\ .$

This leads to the following differential equation for $x(t)$:

$A(t)\dot{x}(t)=-\dot{A}(t)x(t)+\dot{b}(t)-\lambda(A(t)x(t)-b(t))\ .$ (2)

And we have transformed the time-varying linear equations problem into an
initial value differential equations problem that needs to be solved for
$t>0$. This is where the different dynamical systems methods split their ways.
In Zhang Neural Networks, the continuous time differential equation (2) is
then discretized for $0<t_{j}<t_{end}$ and the ensuing derivatives are
approximated by high order difference quotients, with the one for the unknown
$\dot{x}(t_{j})$ being 1-step ahead and proven convergent, while the others
such as for $\dot{A}(t)$ and $\dot{b}(t)$ can be backward difference formulas.
How to proceed from (2) with solving $A(t)x(t)=b(t)$ via ZNN methods is still
an open problem, especially for large scale sparse or structured time-varying
linear equations since the matrix $A(t_{j})$ encumbers the unknown
$\dot{x}(t_{j})$ on the left hand side and there is no known 1-step ahead
differentiation formula that can be used here.
The general idea that underlies ZNN methods for time-varying problems is to
replace repeated matrix computations by solving linear differential equations
and associated initial value problems for discrete instances $0<t_{j}<t_{end}$
instead.

###  2.2 A Zhang Neural Network approach to find time-varying generalized
matrix inverses $\mathbf{Y(t)}$ for time-varying full rank matrices
$\mathbf{B(t)}$ with $\mathbf{B(t)_{m,n}Y(t)_{n,m}=I_{m}}$

This section is based on joint work with Jian Li, Mingzhi Mao, and Yunong
Zhang, see [23].

Continuous Problem Formulation :

For an $m$ by $n$ real time varying matrix $B(t)$ of full rank $m$ with $m\leq
n$ we form the matrix-valued error function

$E(t)=B(t)-Y^{+}(t)\in\mathbb{R}^{m\times n}$ (3)

where the upper + sign always means ’generalized inverse’. Then we use the
design formula

$\dot{E}(t)=-\lambda E(t)$ (4)

with design parameter $\lambda>0$. Based on [24, Lemma 3] we have

$\dot{Y}^{+}(t)=-Y^{+}(t)\dot{Y}(t)Y^{+}(t).$ (5)

And from (3) and (5) we obtain

$\dot{E}(t)=\dot{B}(t)-\dot{Y}^{+}(t)=\dot{B}(t)+Y^{+}(t)\dot{Y}(t)Y^{+}(t).$
(6)

Combining (4) and (6), we then get

$\dot{B}(t)+Y^{+}(t)\dot{Y}(t)Y^{+}(t)=-\lambda(B(t)-Y^{+}(t))$ (7)

And by right multiplying (7) with $Y(t)$ we have

$(\dot{B}(t)+Y^{+}(t)\dot{Y}(t)Y^{+}(t))Y(t)=-\lambda(B(t)-Y^{+}(t))Y(t).$ (8)

With $m\leq n$, we have $Y_{m\times n}^{+}(t)Y_{n\times m}(t)=I_{m\times m}$
and thus

$\dot{B}(t)Y(t)+Y^{+}(t)\dot{Y}(t)=-\lambda(B(t)Y(t)-I).$ (9)

The solution of a generalized-matrix inverse problem is not unique when $m<n$
and we only need to find a solution that satisfies (9). Consequently the
continuous model can be represented as

$\begin{split}\dot{Y}(t)=-\lambda(Y(t)B(t)Y(t)-Y(t))-Y(t)\dot{B}(t)Y(t)\end{split}$
(10)

which agrees completely with [14, formula (15), p. 317]. Substituting (10)
into (9), we have

$\dot{B}(t)Y(t)+Y^{+}(t)(-\lambda(Y(t)B(t)Y(t)-Y(t))-Y(t)\dot{B}(t)Y(t))=-\lambda(B(t)Y(t)-I),$
(11)

which we rewrite as

$\dot{B}(t)Y(t)+(-\lambda(Y^{+}(t)Y(t)B(t)Y(t)-Y^{+}(t)Y(t))-Y^{+}(t)Y(t)\dot{B}(t)Y(t))=-\lambda(B(t)Y(t)-I).$
(12)

With $Y_{m\times n}^{+}(t)Y_{n\times m}(t)=I_{m\times m}$, we have

$\dot{B}(t)Y(t)+(-\lambda(B(t)Y(t)-I)-\dot{B}(t)Y(t))=-\lambda(B(t)Y(t)-I).$
(13)

Thus model (10) satisfies (9) and its solution solves the time-varying
generalized matrix inverse problem.

Zhang Neural Network Discretization :

Given a sequence of rectangular matrices $B_{j}$ at time instances $t_{j}\leq
t_{k}$ we want to find the discrete time-varying generalized matrix inverse
$Y_{k+1}$ of $B_{k+1}$ on each computational time interval
$[k\tau,(k+1)\tau)\subseteq[0,t_{\textnormal{f}}]$ so that

$B_{k+1}-Y^{+}_{k+1}={0}.$ (14)

Here $B_{k+1}=B(t_{k+1})=B((k+1)\tau)\in\mathbb{R}^{m\times n}$ is a time-
varying full rank equidistant matrix sequence, $m\leq n$ and
$Y_{k+1}\in\mathbb{R}^{n\times m}$ is unknown. $Y_{k+1}$ needs to be computed
in real-time for each time interval $[k\tau,(k+1)\tau)\subseteq[0,t_{end}]$.
Here the matrix operator + denotes the generalized inverse of a matrix and
${0}\in\mathbb{R}^{m\times n}$ is the zero matrix. Besides, $k=0,1,\cdots$
denotes the updating index, $t_{end}$ denotes the task duration and $\tau$
denotes the constant sampling gap of the time-varying matrix sequence
$B_{k+1}$. For $m>n$, the procedure is similar.

Note that we must obtain each $Y_{k+1}$ at or before time $t_{k+1}$ for real-
time calculations while the actual value of $B_{k+1}$ is unknown before
$t_{k+1}$. Thus we cannot obtain the solution by calculating
$Y_{k+1}=B^{+}_{k+1}$. To obtain $Y_{k+1}$ in real-time, we must develop a
model based on the available information from before $t_{k+1}$ such as that in
$B_{j}$, $Y_{j}$ and $Y_{j-1}$ for $j\leq k$ instead of unknown information
such as $B_{k+1}$.

To obtain a discrete-time model that solves the original discrete time-varying
generalized matrix inverse problem (14), we need to discretize the continuous
model (10). First we use the conventional 1-step forward Euler formula

$\dot{x}(t_{k})\approx\frac{x(t_{k+1})-x(t_{k})}{\tau}$ (15)

with truncation error of order $O(\tau)$. Based on formula (15) we approximate

$\dot{Y}(t_{k})=\frac{Y(t_{k+1})-Y(t_{k})}{\tau}\ $ (16)

and use this equation to discretize the continuous model (10) as follows

$\begin{split}Y_{k+1}=-h(Y_{k}B_{k}Y_{k}-Y_{k})-\tau
Y_{k}\dot{B}_{k}Y_{k}+Y_{k}.\end{split}$ (17)

Here $h=\tau\lambda$. In most real-world applications, information of the
first-order time derivatives, i.e., the value of $\dot{B}_{k}$ may not be
explicitly known for the discrete time-varying generalized matrix inverse
problem (14). If this is so, the value of $\dot{B}_{k}$ can be approximated by
a backward finite difference formula. To assure the accuracy and simplicity of
the discretized model, the truncation error of the backward finite difference
formula for $\dot{B}_{k}$ should be near equal to that of the 1-step-ahead
finite difference formula that approximates $\dot{Y}_{k}$. Thus $\dot{B}_{k}$
in (17) should best be approximated by Euler’s backward finite difference
formula

$\dot{b}_{k}\approx\frac{b_{k}-b_{k-1}}{\tau},$ (18)

because the truncation error order $O(\tau)$ of formula (18) equals that of
formula (15). Thus we approximately have

$\dot{B}_{k}=\frac{B_{k}-B_{k-1}}{\tau}.$ (19)

Then we combine equation (17) with equation (19) and the Euler discrete model
becomes

$\begin{split}Y_{k+1}=-h(Y_{k}B_{k}Y_{k}-Y_{k})-Y_{k}(B_{k}-B_{k-1})Y_{k}+Y_{k}.\end{split}$
(20)

Note that the truncation error of the discrete model (20) is of order
$\mathbf{O}(\tau^{2})$ where the symbol $\mathbf{O}(\tau^{2})$ denotes a
matrix in which each entry is of order ${O}(\tau^{2})$. This model uses only
present or past information of $B_{k}$, $B_{k-1}$ and $Y_{k}$ and solves for
$Y_{k+1}$. Thus $Y_{k+1}$ can be calculated during the time interval
$[t_{k},t_{k+1})$ and if $Y_{k+1}$ can be computed quickly enough in real-time
it will be ready when time instant $t_{k+1}$ arrives.

Higher accuracy 1-step ahead formulas exist for discrete models, namely

$\dot{x}(t_{k})\approx\frac{2x(t_{k+1})-3x(t_{k})+2x(t_{k-1})-x(t_{k-2})}{2\tau}$
(21)

and

$\dot{x}(t_{k})\approx\frac{6x(t_{k+1})-3x(t_{k})-2x(t_{k-1})-x(t_{k-2})}{10\tau}.$
(22)

Both have truncation errors of order $O(\tau^{2})$. For simplicity we only
consider formula (21) and call it the 4-IFD formula because 4 instants in time
are used to approximate the first-order derivative of $x(t_{k})$. When we
employ the 4-IFD formula (21) inside the our continuous model (10) we obtain

$\begin{split}Y_{k+1}=-h(Y_{k}B_{k}Y_{k}-Y_{k})-\tau
Y_{k}\dot{B}_{k}Y_{k}+\frac{3}{2}Y_{k}-Y_{k-1}+\frac{1}{2}Y_{k-2}.\end{split}$
(23)

Next we use the three-instant backward finite difference formula

$\dot{b}_{k}\approx\frac{3b_{k}-4b_{k-1}+b_{k-2}}{2\tau}$ (24)

with error order $O(\tau^{2})$ to approximate the value of $\dot{B}_{k}$ in
equation (23). Then the 4-IFD-type discretized model becomes

$\begin{split}Y_{k+1}=-h(Y_{k}B_{k}Y_{k}-Y_{k})-Y_{k}\left(\frac{3}{2}B_{k}-2B_{k-1}+\frac{1}{2}B_{k-2}\right)Y_{k}+\frac{3}{2}Y_{k}-Y_{k-1}+\frac{1}{2}Y_{k-2}.\end{split}$
(25)

Its truncation error is of order $\mathbf{O}(\tau^{3})$. Similar to the Euler
based discrete model (20), the 4-IFD-type discrete model uses only present and
past information such as $B_{k}$, $B_{k-1}$, $B_{k-2}$, $Y_{k}$, $Y_{k-1}$ and
$Y_{k-2}$ to solve for $Y_{k+1}$. Thus it also satisfies the requirements for
real-time computation.

A 5 instants finite difference formula :

Any usable finite difference formula for discretizing the continuous model
(10) must satisfy several restrictions. It must be 1-step ahead for $\dot{x}$,
i.e., approximate $\dot{x}(t_{k})$ by using $x_{k+1},\ x_{k},\ x_{k-1}$ and
earlier data, and it must be $0$-stable and convergent. However, 1-step ahead
finite difference formulas do not necessarily generate stable and convergent
discrete models, see e.g., [30, 29].

Here is a new 1-step ahead finite difference formula with higher accuracy than
the Euler and 4-IFD formulas. It will be used to generate a stable and
convergent discrete model that finds time-varying generalized matrix inverses
more accurately in real-time.
Theorem 1  _The 5-IFD formula_

$\dot{x}(t_{k})\approx\frac{8x(t_{k+1})+x(t_{k})-6x(t_{k-1})-5x(t_{k-2})+2x(t_{k-3})}{18\tau}$
(26)

_has truncation error order $O(\tau^{3})$._
The proof relies on four Taylor expansions that use $x(t_{k+1})$ and
$x(t_{k-1})$ through $x(t_{k-3})$ around $x(t_{k})$ and clever linear
combinations thereof.
The new 1-step ahead discretization formula (26) then leads to the 5 instants
discrete model

$\begin{split}Y_{k+1}=&-\frac{9}{4}h(Y_{k}B_{k}Y_{k}-Y_{k})-\frac{9}{4}Y_{k}\left(\frac{11}{6}B_{k}-3B_{k-1}+\frac{3}{2}B_{k-2}-\frac{1}{3}B_{k-3}\right)Y_{k}\\\
&\qquad-\frac{1}{8}Y_{k}+\frac{3}{4}Y_{k-1}+\frac{5}{8}Y_{k-2}-\frac{1}{4}Y_{k-3}.\end{split}$
(27)

which has a truncation error of order $\mathbf{O}(\tau^{4})$.
Theorem 2  _The 5 instants discrete model ( 27) is 0-stable._
The multistep formula of the 5 instants discrete model time-varying
generalized matrix inverses has the characteristic polynomial

$P_{4}(\theta)=\theta^{4}+\frac{1}{8}\theta^{3}-\frac{3}{4}\theta^{2}-\frac{5}{8}\theta+\frac{1}{4}$
(28)

with four distinct roots $\theta_{1,2}=-0.7160\pm 0.5495\textnormal{i}$,
$\theta_{3}=0.3069$ and $\theta_{4}=1$ inside the complex unit circle, making
this model 0-stable.

Numerical Examples :

Example 1. Consider the discrete time-varying generalized matrix inverse
problem

$B_{k+1}-Y^{+}_{k+1}={0},\leavevmode\nobreak\
\textnormal{with}\leavevmode\nobreak\
B_{k}=\begin{bmatrix}\sin(0.5t_{k})&\cos(0.1t_{k})&-\sin(0.1t_{k})\\\
-\cos(0.1t_{k})&\sin(0.1t_{k})&\cos(0.1t_{k})\end{bmatrix}.$ (29)

(a) With $\tau=0.1$ s (b) With $\tau=0.01$ s (c) With $\tau=0.001$ s

Figure 1: Typical residual errors generated by the 5 instants , the 4 instants
and the Euler formulas with different sampling gaps $\tau$ when solving the
discrete time-varying generalized-matrix-inverse problem (29) for $t_{end}=30$
s, and $h=0.1$.

Example 2. Here we consider the discrete time-varying matrix inversion problem

$A_{k+1}X_{k+1}=I\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \textnormal{with}\leavevmode\nobreak\
A_{k}=\begin{bmatrix}\sin(0.5t_{k})+2&\cos(0.5t_{k})\\\
\cos(0.5t_{k})&\sin(0.5t_{k})+2\\\ \end{bmatrix}.$ (30)

(a) Profiles of the (1,1) entry $X_{1,1}$ (b) Profiles of the (1,2) entry
$X_{1,2}$ (c) Profiles of the (2,1) entry $X_{2,1}$ (d) Profiles of the (2,2)
entry $X_{2,2}$

Figure 2: Profiles of the four entries of the solution $X$ when solving the
discrete time-varying matrix-inverse problem (30) with $\tau=0.1$ s. Here the
solid curves show the solution entries generated by the 5 instants discrete
model obtained from random starting values and the dash-dotted curves depict
the theoretical solutions.

###  2.3 A Zhang Neural Network approach for solving nonlinear convex
optimization problems under time-varying linear constraints

This section is based on joint work with Jian Li, Mingzhi Mao, and Yunong
Zhang, see [22].

Problem formulation :

$\begin{split}&\text{find }\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{ min}\leavevmode\nobreak\ \leavevmode\nobreak\
f(x(t),t)\\\ &\text{such that }\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ A(t)x(t)=b(t)\
\ \ \text{ with }x(t)\in\mathbb{R}^{n},\ b(t)\in\mathbb{R}^{m}\text{ and
}A(t)\in\mathbb{R}^{m,n}.\end{split}$

Building a continuous time model, i.e., formulating an equation for the
problem :
The Zhang Neural Network approach can be built on the Lagrange function

$L(x(t),l(t),t)=f(x(t),t)+l^{T}(t)(A(t)x(t)-b(t)),$

where $l(t)=[l_{1}(t),\cdots,l_{m}(t)]^{T}\in\mathbb{R}^{m}$ is the Lagrange-
multiplier vector and $..^{T}$ denotes the transpose. Note that there will no
need to solve for the Lagrange functions $l(t)$ here. Set

$y(t)=[x^{T}(t),l^{T}(t)]^{T}=[y_{1}(t),\cdots,y_{n}(t),y_{n+1}(t),\cdots,y_{n+m}(t)]^{T}\in\mathbb{R}^{n+m}$

and

$\begin{split}h(y(t),t)=\begin{bmatrix}\dfrac{\partial f(x(t),t)}{\partial
x}+A^{T}(t)l(t)\\\\[8.53581pt]
A(t)x(t)-b(t)\end{bmatrix}=\begin{bmatrix}h_{1}(y(t),t)\\\ \vdots\\\
h_{n}(y(t),t)\\\ h_{n+1}(y(t),t)\\\ \vdots\\\
h_{n+m}(y(t),t)\end{bmatrix}\in\mathbb{R}^{n+m}.\end{split}$

We transform the multiplier problem into an initial value DE problem instead.
By stipulating exponential decay for $h(t)$ we obtain the model equation

$\dot{y}(t)=-H^{-1}(y(t),t)\left(\lambda hy(t),t)+\dot{h}_{t}(y(t),t)\right)$

for the Jacobian matrix

$H(y(t),t)=\begin{bmatrix}\dfrac{\partial f^{2}(x(t),t)}{\partial
x}\partial^{T}x&&A^{T}(t)\\\\[8.53581pt] A(t)&&0\end{bmatrix}\ \text{ and }\
\dot{h}_{t}(y(t),t)=\dfrac{\partial h(y(t),t)}{\partial t}\ .$

Discretizing the model and choosing suitable high order finite difference
formulas :
To discretize the continuous model

$\begin{split}\dot{y}(t)=-H^{-1}(y(t),t)\left(\lambda
hy(t),t)+\dot{h}_{t}(y(t),t)\right)\end{split}$

we can use the forward Euler difference formula with truncation error order
$O(\tau)$

$\dot{x}(t_{k})=\frac{x(t_{k+1})-x(t_{k})}{\tau}$

or the four-instance forward difference formula (4-IFD)

$\dot{x}(t_{k})=\frac{5x(t_{k+1})-3x(t_{k})-x(t_{k-1})-x(t_{k-2})}{8\tau}$

with truncation error order $O(\tau^{2})$. The Euler formula yields the
discretized model

$y_{k+1}=-H^{-1}(y_{k},t_{k})\left(\kappa
h(y_{k},t_{k})+\tau\dot{h}_{t}(y_{k},t_{k})\right)+y_{k}\ \text{ with }\
\kappa=\tau\lambda$

while the 4-IFD formula results in

$y_{k+1}=-\frac{8}{5}H^{-1}(y_{k},t_{k})\left(\kappa
h(y_{k},t_{k})+\tau\dot{h}_{t}(y_{k},t_{k})\right)+\frac{3}{5}y_{k}+\frac{1}{5}y_{k-1}+\frac{1}{5}y_{k-2}+O(\tau^{3})\
.$

Both discretization formulas are consistent and convergent. This can be proved
via the the roots of the associated characteristic polynomial. Its roots must
lie in the complex unit circle and cannot be repeated on its boundary.
Since the value of $\dot{h}_{t}(y_{k},t_{k})$ may not be known explicitly we
may replace it by

$\dot{h}_{t}(y_{k},t_{k})=\frac{3h(y_{k},t_{k})-4h(y_{k},t_{k-1})+h(y_{k},t_{k-2})}{2\tau}$

which uses the three-point backward finite difference formula

$\dot{x}(t_{k})=\frac{3x(t_{k})-4x(t_{k-1})+x(t_{k-2})}{2\tau}\ $

of order $O(\tau^{2})$.
Then the discretized FIFD formula becomes more complicated but easier to
implement:

$\displaystyle y_{k+1}$ $\displaystyle=$
$\displaystyle-\frac{8}{5}H^{-1}(y_{k},t_{k})\left(\left(\kappa+\frac{3}{2}\right)h(y_{k},t_{k})-2h(y_{k},t_{k-1})+\frac{1}{2}h(y_{k},t_{k-2})\right)$
$\displaystyle\ \ \ +\frac{3}{5}y_{k}+\frac{1}{5}y_{k-1}+\frac{1}{5}y_{k-2}\
\quad\text{ of order }\ \ O(\tau^{3})\ .$

To implement this formula, the inverse of the Jacobian matrix $H$ can be
computed at each time $t_{k}$ in a fraction of the available real-time
interval $[t_{k},t_{k+1})$ by using the real-time inverse finding ZNN method
from the previous subsection 2.2.

Numerical example and results :
As an example we solve the following convex nonlinear optimization problem
with known theoretical solution numerically by using our ZNN method, for
further details and applications see [22]:

$\begin{split}\text{Find }\ \
&\text{min}(\cos(0.1t_{k+1})+2)x_{1}^{2}+(\cos(0.1t_{k+1})+2)x_{2}^{2}+2\sin(t_{k+1})x_{1}x_{2}+\sin(t_{k+1})x_{1}+\cos(t_{k+1})x_{2}\\\
&\text{so that}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\sin(0.2t_{k+1})x_{1}+\cos(0.2t_{k+1})x_{2}=\cos(t_{k+1}).\end{split}$

The FIFD formula is a 4 instance formula, while the Euler formula needs only
two. Both discretization models work in real-time and both typically create
the optimal solution in a fraction of a second with differing degrees of
accuracy according to their orders.
The example below runs for 10 sec. The time-varying values for $f(x(t),t)$,
$A(t)$ and $b(t)$ are given as functions and evaluated from their function
formulations. In real world applications these values might be supplied by
sensors during each time interval $t_{i}\leq t_{i+1}$ and the empirical values
would be inserted into the difference formulas as they are evaluated by
sensors in real time.

Figure 3:  Solution states with $\tau=0.1$ sec, and solution errors generated
by the 4-IFD based discretization model (of order $O(\tau^{3})$) and the Euler
based model (of order $O(\tau^{2})$) with $\tau=0.1$, $0.01$, and $0.001$ sec
and $\lambda=10$:
(a) solution states with $\tau=0.1$ sec, (b) solution errors with $\tau=0.1$
sec,
(c) solution errors with $\tau=0.01$ sec, (d) solution errors with
$\tau=0.001$ sec.

##  3 On Quantum and Multi-state Computing (_Epoch 8_ yet to start and come)

Quantum Computing and multi-state memory and computers with multi-state
processors will change the way we compute once they become available. They
will require new operating systems and new software with new and yet to be
discovered algorithms. What will this new era entail? Nobody knows or can
reliably predict.
I asked an ‘expert’ on quantum computing three years ago as to when he
expected to have a quantum computer at his disposal or on his desk. The answer
was : “Not in my lifetime, not in 20 years”.
Currently about a dozen or more research centers in Europe and South-East Asia
are trying to build quantum computers based on the quantum superposition
principle and quantum entanglement of elementary particles. They do so in a
multitude of different ways. The envisioned benefit of these efforts would be
to be able to compute super fast in parallel and in simulations to solve huge
data problems quicker than ever before and to solve problems that are
unassailable now with our current best supercomputer networks. All of the
proposed quantum science techniques make use of superconducting circuits and
particles. The aim is to build quantum computers in one or two decades with
around 100 entangled quantum bits. Such a quantum computer would be bulky, it
would need much supplementary equipment for cooling and so forth and could
easily take up a whole floor of a building, just as the first German and
British valve computers did in the 1940s. But it would surpass the computing
capacity of all current supercomputers and desk and laptops on Earth combined.
Currently the largest working entangled quantum array contains fewer than 10
quantum bits. Access of a 100 bit quantum computer would probably be via the
Cloud and there would be no quantum computer laptops. Quantum computers may
take another 10, 20 or 30 years to materialize.
_How will they come about? Which yet unknown algorithms will they use? Who
will invent them? Who code them?_
If history can be a guide, John Francis and Vera Kublanovskaya were both
working independently on circuit diagrams and logic gate designs for valve
computers in England and in Russia at the time when they discovered QR (or LQ)
in the late 1950s.
So we possibly are looking for quantum computer hardware and software
designers who know numerical analysis and algorithm development in or about
the year 2040. In a similar fashion, Leibnitz and Seki formalized our now
ubiquitous matrix concept independently but simultaneously in 1683, in Germany
and in Japan.
_Maybe it will take two again?_

The references given below only go back to the year 1799.

## References

  * [1]
  * [2] N. H. Abel, Mémoire sur les équations algébriques où on démontre l’impossibilité de la résolution de l’équation générale du cinquième degré, Groendahl, Christiana (Kopenhagen) (1824), 7 p.
  * [3] W. E. Arnoldi, The principle of minimized iterations in the solution of the matrix eigenvalue problem, Quarterly of Applied Mathematics, volume 9 (1951), p. 17 - 29.
  * [4] K. Braman, R. Byers, R. Mathias, The Multishift QR Algorithm. Part I: Maintaining Well-Focused Shifts and Level 3 Performance, SIAM J. Matrix Anal. & Appl., 23 (2006), 929 - 947.
  * [5] K. Braman, R. Byers, R. Mathias, The Multishift QR Algorithm. Part II: Aggressive Early Deflation, SIAM J. Matrix Anal. & Appl., 23 (2006), 948 - 973.
  * [6] A. L. Cauchy, Sur l’equation à l’aide de laquelle on determine les inégalités séculaires des mouvements des planètes, Exerc. de Math., 4, (1829); also Oeuvres (2) 9, p. 174-195.
  * [7] A. L. Cauchy, Mémoire sur l’intégration des équations linéaires, Comptes Rendus, 8, (1839), p. 827-830, 845-865, 889-907, 931-937.
  * [8] J. G. F. Francis, The QR Transformation, I, The Computer Journal, 4 (1961), p. 265 - 271.
  * [9] J. G. F. Francis, The QR Transformation, II, The Computer Journal, 4 (1962), p. 332 - 345.
  * [10] É. Galois, Analyse d’un mémoire sur la résolution algébraique des équations, Bulletin des Sciences mathématiques, physiques et chimiques, XIII (1830), p. 271-272.
  * [11] É. Galois, Note sur la résolution des équations numériques, Bulletin des Sciences mathématiques, physiques et chimiques, XIII (1830), p. 413-414.
  * [12] É. Galois, Mémoire sur les conditions de résolubilité des équations par radicaux. Journal de mathématiques pures et appliquées, (published by Joseph Liouville), XI (1846), p. 417-433.
  * [13] C. F. Gauß , Demonstratio nova theorematis omnem functionem algebraicam rationalem integram unius variabilis in factores reales primi vel secundi gradus resolvi posse, PhD thesis, Universität Helmstedt, (1799), Werke III, p. 1 30.
  * [14] Neil Getz and ]errold E. Marsden, Dynamical methods for polar decomposition and inversion of matrices, Linear Algebra and its Appl., 258 (1997), p. 311 - 343.
  * [15] Gene Golub, W. Kahan, Calculating the singular values and pseudo-inverse of a matrix, J. SIAM Numer. Anal. Ser. B, 2 (1965), p. 205 - 224.
  * [16] G. Golub, F. Uhlig, The QR algorithm: 50 years later; its genesis by John Francis and Vera Kublanovskaya, and subsequent developments, IMA J. Num. Anal., vol. 29 (2009), p. 467 - 485.
  * [17] W. G. Horner, A new method of solving numerical equations of all orders, by continuous approximation, Philosophical Transactions. Royal Society of London, (July 1819), p. 308 - 335.
  * [18] C. G. Jacobi, (1839), Canon Arithmeticus, Berlin, Typis Academicis Berolini, 1839.
  * [19] V. N. Kublanovskaya, On some algorithms for the solution of the complete eigenvalue problem, USSR Computational Mathematics and Mathematical Physics, 1 (1963), p. 637 - 657 (1963, received Feb 1961). Also published in: Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki [Journal of Computational Mathematics and Mathematical Physics], 1 (1961), p. 555 - 570.
  * [20] A. N. Krylov, On the Numerical Solution of Equation by Which are Determined in Technical Problems the Frequencies of Small Vibrations of Material Systems, News of Academy of Sciences of the USSR), Otdel. mat. i estest. nauk, VII, Nr.4 (1931), p. 491 - 539 (in Russian).
  * [21] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Res. Nat. Bur. Stand., 45 (1950), p. 255 - 282.
  * [22] Jian Li, Mingzhi Mao, Frank Uhlig, Yunong Zhang, Z-type neural-dynamics for time-varying nonlinear optimization under a linear equality constraint with robot application, Journal of Computational and Applied Mathematics, 327 (2018), p. 155-166.
  * [23] Jian Li, Mingzhi Mao, Frank Uhlig, Yunong Zhang, A 5-Instant finite difference formula to find discrete time-varying generalized matrix inverses, matrix inverses and scalar reciprocals, 17 p, submitted.
  * [24] B. Liao, Y. Zhang, From different ZFs to different ZNN models accelerated via Li activation functions to finite-time convergence for time-varying matrix pseudoinversion, Neurocomputing 133 (2014), p. 512 - 522.
  * [25] Heinz Rutishauser, Une m ethode pour la d etermination des valeurs propres d une matrice, C. R. Acad. Sci., 240 (1955), p. 34 - 36.
  * [26] Heinz Rutishauser, Solution of eigenvalue problems with the LR-transformation, Nat. Bur. Stand. Appl. Math. Ser., 49 (1958), p. 47 - 81.
  * [27] Olga Taussky, Note on the condition of matrices, Math. Tables and Other Aids to Computation, 4 (1950), p. 111 - 112.
  * [28] Yunong Zhang, Danchi Jiang, Jun Wang, A recurrent neural network for solving Sylvester equation with time-varying coefficients, IEEE Transactions on Neural Networks 13 (2002), p. 1053 - 1063.
  * [29] Y. Zhang, Y. Chou, J. Chen, Z. Zhang, and L. Xiao, Presentation, error analysis and numerical experiments on a group of 1-step-ahead numerical differentiation formulas, J. Comput. Appl. Math., 239 (2013), pp. 406–414.
  * [30] Y. Zhang, L. Jin, D. Guo, Y. Yin, and Y. Chou, Taylor-type 1-step-ahead numerical differentiation rule for first-order derivative approximation and ZNN discretization, J. Comput. Appl. Math., 273 (2014), p. 29 - 40.