# High-order parallel-in-time method for the monodomain equation in cardiac
electrophysiology††thanks: This work was supported by the European High-
Performance Computing Joint Undertaking EuroHPC under grant agreements No
955495 (MICROCARD) and No 955701 (TIME-X), co-funded by the Horizon 2020
programme of the European Union (EU) and the Swiss State Secretariat for
Education, Research and Innovation.

Giacomo Rosilho de Souza Corresponding author<EMAIL_ADDRESS>Euler Institute, Università della Svizzera italiana, via G. Buffi 13, Lugano,
6900, Switzerland  Simone Pezzuto Contributing author<EMAIL_ADDRESS>Laboratory for Mathematics in Biology and Medicine, Department of Mathematics,
Università di Trento, via Sommarive 14, Trento, 38123, Italy  Rolf Krause
Contributing author<EMAIL_ADDRESS>Faculty of Mathematics and
Informatics, FernUni, Schinerstrasse 18, Brig, 3900, Switzerland  AMCS, CSE
Department, KAUST, Thuwal-Jeddah, Saudi Arabia

###### Abstract

Simulation of the monodomain equation, crucial for modeling the heart’s
electrical activity, faces scalability limits when traditional numerical
methods only parallelize in space. To optimize the use of large multi-
processor computers by distributing the computational load more effectively,
time parallelization is essential. We introduce a high-order parallel-in-time
method addressing the substantial computational challenges posed by the stiff,
multiscale, and nonlinear nature of cardiac dynamics. Our method combines the
semi-implicit and exponential spectral deferred correction methods, yielding a
hybrid method that is extended to parallel-in-time employing the PFASST
framework. We thoroughly evaluate the stability, accuracy, and robustness of
the proposed parallel-in-time method through extensive numerical experiments,
using practical ionic models such as the ten-Tusscher–Panfilov. The results
underscore the method’s potential to significantly enhance real-time and high-
fidelity simulations in biomedical research and clinical applications.

Keywords— cardiac electrophysiology, monodomain equation, parallel-in-time,
PFASST, spectral deferred correction.

MSC Classification— 65L04,65L10,65L20,65Y05.

##  1 Introduction

The heart’s ability to pump blood hinges on the propagation of electrical
signals across cardiac cells and tissues. Central to modeling these complex
dynamics is the monodomain equation, which characterizes the propagation of
the electric potential and is modulated by ionic models. The latter, describe
the interaction of ionic currents and cellular membranes, encompassing a broad
range of physiological and pathological cardiac behaviors.

Simulating the monodomain equation is computationally demanding due to its
multiscale, nonlinear nature and the high spatio-temporal resolutions required
for capturing cardiac dynamics. In traditional numerical methods, optimized
for spatial parallel scalability, time integration can be a bottleneck.
Indeed, to increase accuracy, both spatial and time resolutions are improved.
Although increasing the number of spatial processors keeps the cost per time
step constant, it doesn’t address the need for additional time steps, and the
number of processors is limited by communication bottlenecks. As a result,
these schemes have limited scalability, impeding real-time simulations. To
mitigate these limitations, we introduce a high-order parallel-in-time (PinT)
method specifically for the monodomain equation. This method enhances
scalability by allowing parallelization in the time direction, especially
beneficial when spatial scalability is maxed out.

The literature on PinT methods for the monodomain equation is very scarce. In
[1], the authors solve the monodomain equation in parallel-in-time employing
the parareal framework [28]. In [2], the PFASST [15] and Space-Time Multigrid
[19] methods are compared on a reaction-diffusion problem, which turns out to
be a special case of the monodomain equation. In both cases the scalar
nonstiff FitzHugh–Nagumo (FHN) [17, 30] ionic model is employed, which is
overly simplistic to be employed in practical simulations. Due to stability
and efficiency reasons, the methods in [1, 2] cannot be used with realistic
highly nonlinear and multi-scale models, such as the popular ten-
Tusscher–Panfilov (TTP) [41] or the Courtemance–Ramirez–Nattel (CRN) [12]
models, for instance.

After spatial discretization, the monodomain equation can be written as

$y^{\prime}=f_{I}(t,y)+f_{E}(t,y)+f_{e}(t,y),\qquad y(0)=y_{0},$ (2a) where
$f_{e}(t,y)=\Lambda(y)(y-y_{\infty}(y)).$ (2b)

In 2, $f_{I}$ is a stiff linear term representing diffusion of the electric
potential, $f_{e}$ is very stiff and drives the dynamics of the ionic model’s
gating variables, and $f_{E}$ is a nonstiff but nonlinear and expensive to
evaluate term for the remaining ionic model’s variables. The stiff matrix
$\Lambda(y)$ is diagonal and $y_{\infty}(y)$ is a nonlinear term for the
gating variables’ equilibrium.

The method proposed here leverages the PFASST framework, which integrates the
parareal [28], spectral deferred correction (SDC) [14], and full approximation
scheme (FAS) [42] methods. The parareal algorithm, an iterative PinT method,
starts with a coarse solver that progresses through the entire time domain
serially; quickly but inaccurately. Then, in parallel, fine solvers refine the
results in smaller sub-intervals. Enhancements to parareal incorporate the SDC
method for both fine and coarse solvers. SDC is an iterative method for
solving the monolithic nonlinear system in implicit Runge–Kutta (RK) methods.
It works as a preconditioned fixed point iteration and is flexible: able to
adjust its preconditioner for specific problems, apply to any RK method, and
achieve arbitrary accuracy. When employed in the parareal framework, SDC
methods of high and low order serve as fine and coarse solvers, respectively.
Furthermore, SDC synergizes with parareal by blending their iterations and
improving parallel efficiency. Combined with FAS to ensure precise information
transfer between levels, this gives rise to the multi-level PinT PFASST method
[15, 16, 3].

Based on this framework, new PinT methods are designed by replacing SDC with
one of its variants. For instance, the semi-implicit SDC (SI-SDC) [29], for
problems 2 with $f_{e}\equiv 0$. Or the exponential SDC (ESDC) [6], for 2 with
$f_{I}\equiv f_{E}\equiv 0$. In this paper, we propose a new hybrid SDC (HSDC)
method that combines SI-SDC, for $f_{I}$ and $f_{E}$, with ESDC, for $f_{e}$.
A PinT method based on HSDC is seamlessly obtained using the PFASST
methodology. For simplicity, we will denote the method HSDC in both its serial
and parallel versions.

Although the HSDC method is designed for the monodomain equation, it can be
applied to any problem structured as 2. Like SI-SDC and ESDC, it has an
arbitrary order of convergence. This is uncommon for the monodomain equation
even in the serial context; indeed, for efficiency and stability, most methods
employ splitting strategies and are low order. Further, due to exponential
integration, the method is stable even for severely stiff ionic models, such
as TTP. Finally, the method is computationally cost-effective, as its
preconditioner relies on the efficient Rush–Larsen method [32, 27] popular in
cardiac electrophysiology; that is, a splitting method based on implicit,
explicit, and exponential Euler steps. The HSDC method proposed here has been
implemented in pySDC [35], a Python package for prototyping and test SDC-like
methods.

The paper is structured as follows. Section 2 reviews the monodomain equation
and Section 3 describes the development of the PinT HSDC method. In Section 4
we perform a numerical stability analysis of the HSDC method. In Section 5, we
present a comprehensive set of numerical experiments that demonstrate the
robustness and effectiveness of HSDC. Conclusions are drawn in Section 6.

##  2 The monodomain equation for cardiac electrophysiology

The monodomain equation is a fundamental tool in cardiac electrophysiology,
describing how electrical signals propagate through heart tissue. It models
the interplay of ion channels, cellular membranes, tissue conductivity and
electric potential propagation. Therefore, it provides insights into phenomena
ranging from normal cardiac rhythm to pathological conditions such as
arrhythmia.

###  2.1 The monodomain model

The monodomain equation is given by the following reaction-diffusion system of
equations

$\displaystyle C_{m}\frac{\partial V_{m}}{\partial t}$
$\displaystyle=\chi^{-1}\nabla\cdot\bigl{(}\sigma\nabla
V_{m}\bigr{)}+I_{\text{stim}}(t,\bm{x})-I_{\text{ion}}(V_{m},{w}_{a},{w}_{g}),\qquad$
in $\Omega\times(0,T]$, (3a) $\displaystyle\frac{\partial{w}_{a}}{\partial t}$
$\displaystyle=h_{a}(V_{m},{w}_{a},{w}_{g}),$ in $\Omega\times(0,T]$, (3b)
$\displaystyle\frac{\partial{w}_{g}}{\partial t}$
$\displaystyle=h_{g}(V_{m},{w}_{g})={\alpha}(V_{m})\bigl{(}\mathbf{1}-{w}_{g}\bigr{)}-{\beta}(V_{m}){w}_{g},$
in $\Omega\times(0,T]$, (3c) $\displaystyle-\sigma\nabla V_{m}\cdot\bm{n}$
$\displaystyle=0,$ $\displaystyle\text{on }\partial\Omega\times(0,T],$ (3d)
$\displaystyle V_{m}(0,\bm{x})$
$\displaystyle=V_{m,0},\quad{w}_{a}(0,\bm{x})={w}_{a,0},\quad{w}_{g}(0,\bm{x})={w}_{g,0},$
in $\Omega$ (3e)

where $T>0$ and $\Omega\subset\mathbb{R}^{d}$ is a domain. The model unknowns
are the the transmembrane potential $V_{m}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{d}\times\mathbb{R}\rightarrow\mathbb{R}$, in
$\mathrm{mV}$, the vector of auxiliary ionic variables
$w_{a}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{d}\times\mathbb{R}\rightarrow\mathbb{R}^{m_{1}}$, and
the vector of ionic gating variables $w_{g}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{d}\times\mathbb{R}\rightarrow\mathbb{R}^{m_{2}}$. The
values for the isotropic conductivity coefficient $\sigma$, the surface area-
to-volume ratio $\chi$ and the cellular membrane electrical capacitance
$C_{m}$ are given in Section 5.

The $I_{\text{stim}}$ term is used to stimulate the tissue and initiate an
action potential’s propagation. The ionic current $I_{\text{ion}}$ models the
transmembrane currents and, like $h_{a}$, $\alpha$, $\beta$ and $m_{1}$,
$m_{2}$, its exact definition depends on the chosen ionic model. Several
models exists in the literature, and the choice depends on the modeled
species, type of muscle and pathological conditions, for instance.
Nevertheless, in most cases the structure is as in 3, hence with
$\mathbf{1}\in\mathbb{R}^{m_{2}}$ a vector of ones, $\alpha(V_{m})$,
$\beta(V_{m})$ diagonal matrices and $h_{a}$, $I_{\text{ion}}$ general
nonlinear terms. Therefore,

$\frac{\partial w_{g}}{\partial
t}=\Lambda_{g}(V_{m})\bigl{(}{w}_{g}-{w}_{g,\infty}(V_{m})\bigr{)},$ (4)

where ${\Lambda}_{g}(V_{m})=-({\alpha}(V_{m})+{\beta}(V_{m}))$ is a diagonal
matrix and
${w}_{g,\infty}(V_{m})=({\alpha}(V_{m})+{\beta}(V_{m}))^{-1}{\alpha}(V_{m})\mathbf{1}$
is the steady state value of ${w}_{g}$ (for fixed $V_{m}$).

###  2.2 Semi-discrete equation

For the purpose of this paper we consider rectangular domains $\Omega$ and
discretize the monodomain equation 3 using a fourth-order finite difference
scheme, yielding the semi-discrete system

$\displaystyle
C_{m}\frac{\operatorname{d\\!}\bm{V}_{m}}{\operatorname{d\\!}t}$
$\displaystyle=\chi^{-1}\mathbf{A}\bm{V}_{m}+\bm{I}_{\text{stim}}(t)-\bm{I}_{\text{ion}}(\bm{V}_{m},\bm{w}_{a},\bm{w}_{g}),$
(5a) $\displaystyle\frac{\operatorname{d\\!}\bm{w}_{a}}{\operatorname{d\\!}t}$
$\displaystyle=\bm{h}_{a}(\bm{V}_{m},\bm{w}_{a},\bm{w}_{g}),$ (5b)
$\displaystyle\frac{\operatorname{d\\!}\bm{w}_{g}}{\operatorname{d\\!}t}$
$\displaystyle=\bm{\Lambda}_{g}(\bm{V}_{m})(\bm{w}_{g}-\bm{w}_{g,\infty}(\bm{V}_{m})),$
(5c)

plus the initial conditions, where $\bm{V}_{m}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\rightarrow\mathbb{R}^{N_{dof}}$,
$\bm{w}_{a}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\rightarrow\mathbb{R}^{m_{1}N_{dof}}$,
$\bm{w}_{g}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\rightarrow\mathbb{R}^{m_{2}N_{dof}}$ and $N_{dof}$ is the
number of mesh degrees of freedom. The ionic model’s terms
$\bm{h}_{a},\bm{\Lambda}_{g},\bm{w}_{g,\infty},\bm{I}_{\text{ion}}$ and the
stimulus $\bm{I}_{\text{stim}}$ are evaluated node-wise.

The semi-discrete equation 5 can be written as 2 for
$y=(\bm{V}_{m},\bm{w}_{a},\bm{w}_{g})\in\mathbb{R}^{N}$, with
$N=(1+m_{1}+m_{2})N_{dof}$, by defining

$\displaystyle f_{I}(t,y)=\begin{pmatrix}(\chi
C_{m})^{-1}\mathbf{A}\bm{V}_{m}\\\ \bm{0}\\\ \bm{0}\end{pmatrix},\qquad
f_{E}(t,y)=\begin{pmatrix}C_{m}^{-1}(\bm{I}_{stim}(t)-\bm{I}_{ion}(\bm{V}_{m},\bm{w}_{a},\bm{w}_{g}))\\\
\bm{h}_{a}(\bm{V}_{m},\bm{w}_{a},\bm{w}_{g})\\\ \bm{0}\end{pmatrix},$ (6a) and
$f_{e}(t,y)=\Lambda(y)(y-y_{\infty}(y))$, where
$\Lambda(y)=\begin{pmatrix}\bm{0}&\bm{0}&\bm{0}\\\ \bm{0}&\bm{0}&\bm{0}\\\
\bm{0}&\bm{0}&\bm{\Lambda}_{g}(\bm{V}_{m})\end{pmatrix},\quad
y_{\infty}(y)=\begin{pmatrix}\bm{0}\\\ \bm{0}\\\
\bm{w}_{g,\infty}(\bm{V}_{m})\end{pmatrix}$ (6b)

and $\bm{\Lambda}_{g}(\bm{V}_{m})$ is a diagonal matrix.

Concerning the numerical integration of 5 written in the form of 2, note that
$f_{I}$ is stiff but linear and $\mathbf{A}$ is symmetric, hence it can be
integrated relatively efficiently by an implicit method. The $f_{E}$ term
houses the most intricate components of the ionic model, hence it is
expensive. Yet, it is nonstiff and can therefore be integrated explicitly.
Finally, $f_{e}$ is very stiff but due to the diagonal form of $\Lambda(y)$ it
is easily integrated exponentially. These considerations motivate the
Rush–Larsen method, proposed in [32], where a splitting strategy is used and
$f_{I}$, $f_{E}$, and $f_{e}$ are integrated by the implicit, explicit and
exponential Euler methods, respectively. This approach is still very popular
nowadays, due to its robustness, efficiency and simple implementation [27]. In
Section 3.1.3 we propose a preconditioner based on this approach.

##  3 Parallel-in-time method for the monodomain equation

In this section, we present the PinT method for the monodomain equation
defined by 2 and 6. To design such method, we first formulate a new SDC method
where the three terms $f_{I}$, $f_{E}$, and $f_{e}$ are integrated with
appropriate schemes. Then we derive its PinT extension using the PFASST
framework, to do so we follow the established methodology [4, 3, 33, 34, 35]
employing preconditioned Picard iterations and the multigrid perspective.

###  3.1 Hybrid spectral deferred correction method

Here we derive the HSDC method for the monodomain equation. This method is an
hybrid between ESDC [6] and SI-SDC [29], where the exponential term $f_{e}$ is
integrated with ESDC and the implicit and explicit terms $f_{I},f_{E}$ are
integrated with the SI-SDC method. The construction of such method follows
three basic steps. First, the monodomain equation is written in a form more
compatible with exponential Runge–Kutta (ERK) methods [24, 23]. Then ERK
methods of collocation type [23] are re-derived, we do that for completeness
and mostly to obtain a formulation closer to collocation RK methods and
therefore easily implementable in available software for PinT methods, such as
pySDC [35]. Finally, we recall the preconditioned Picard iteration for solving
the nonlinear system arising in each step of the ERK method and propose a
preconditioner tailored to the monodomain equation.

####  3.1.1 Monodomain equation reformulation

Here we rewrite the monodomain equation in a ERK-friendly formulation.

Let $\Delta t>0$ be the step size and $t_{n}=n\Delta t$. In the time interval
$t\in[t_{n},t_{n+1}]$, equation 2 with initial value $y_{n}$ can be written as

$y^{\prime}=\Lambda(y_{n})(y-y_{n})+g(y),\qquad y(t_{n})=y_{n},$ (7)

where

$g(y)=f_{I}(t,y)+f_{E}(t,y)+f_{e}(t,y)+\Lambda(y_{n})(y_{n}-y)$ (8)

and, for simplicity, $y_{n}$ and the time variable $t$ are dropped from the
arguments of $g$.

In the next section we define the ERK methods for equation 2 written in the
form 7, where the matrix exponential of $\Lambda(y_{n})$ will be needed.
Concerning accuracy, $\Lambda(y_{n})$ could be replaced by any matrix;
however, for stability reasons, this matrix should capture the stiffness of
the problem. In our case we have chosen $\Lambda(y_{n})$, so that it captures
the stiffness of the exponential term $f_{e}$ evaluated at values $y$
sufficiently close to $y_{n}$. In contrast, the stiffness of the $f_{I}$ term
will be resolved implicitly. This choice is particularly efficient since due
to the diagonal form of $\Lambda(y_{n})$, leading to cheap matrix exponential.
Additionally, $f_{I}$ is linear, hence its implicit integration is as well
relatively inexpensive.

####  3.1.2 Exponential Runge–Kutta methods

Here we derive ERK methods of collocation type for 7, but choose a formulation
resembling collocation RK methods.

Applying the variation of constants formula to 7 we obtain111Notice that,
generally, there is an exponential in front of $y_{n}$ [24, 23]. Here, we
prefer to put the exponential inside the integral and take out the identity
operator. This approach leads to a formulation of ERK closer to standard RK
methods.

$y(t)=y_{n}+\int_{t_{n}}^{t}e^{(t-s)\Lambda(y_{n})}g(y(s))\operatorname{d\\!}s,$
(9)

and performing the change of variables $t=t_{n}+r\Delta t$ yields

$y(t_{n}+r\Delta t)=y_{n}+\Delta t\int_{0}^{r}e^{(r-s)\Delta
t\Lambda(y_{n})}g(y(t_{n}+s\Delta t))\operatorname{d\\!}s.$ (10)

To discretize 10 and obtain exponential integrators, $g(y(t_{n}+s\Delta t))$
is approximated by an interpolating polynomial. Let $0<c_{1}<\ldots<c_{M}=1$
be Radau IIA collocation nodes and, for notational purposes, we also define
$c_{0}=0$. We denote
$\bm{y}_{n}=(y_{n,0},\ldots,y_{n,M})\in\mathbb{R}^{N(M+1)}$, where the
$y_{n,j}\in\mathbb{R}^{N}$ represent approximations to $y(t_{n,j})$, with
$t_{n,j}=t_{n}+c_{j}\Delta t$. Let

$\displaystyle p_{g}(\bm{y}_{n})(s)$
$\displaystyle=\sum_{j=1}^{M}g(y_{n,j})\ell_{j}(s),$ with
$\displaystyle\ell_{j}(s)$ $\displaystyle=\prod_{\begin{subarray}{c}k=1\\\
k\neq j\end{subarray}}^{M}\frac{s-c_{k}}{c_{j}-c_{k}},$ (11)

be a polynomial interpolating $g(y_{n,j})$, $j=1,\ldots,M$. Note that
$y_{n,0}$ and $c_{0}$ are excluded from the definitions of $p_{g}(\bm{y}_{n})$
and the Lagrange basis functions $\ell_{j}$. Replacing $g(y(t_{n}+s\Delta t))$
with $p_{g}(\bm{y}_{n})(s)$ in 10 and evaluating at $r=c_{i}$ yields the
system of equations for $y_{n,i}$:

$y_{n,i}=y_{n}+\Delta t\int_{0}^{c_{i}}e^{(c_{i}-s)\Delta
t\Lambda(y_{n})}p_{g}(\bm{y}_{n})(s)\operatorname{d\\!}s,\qquad i=0,\ldots,M.$
(12)

Let

$a_{ij}(z)=\int_{0}^{c_{i}}e^{(c_{i}-s)z}\ell_{j}(s)\operatorname{d\\!}s,$
(13)

then 12 is rewritten as

$y_{n,i}=y_{n}+\Delta t\sum_{j=1}^{M}a_{ij}(\Delta
t\Lambda(y_{n}))g(y_{n,j}),\qquad i=0,\ldots,M,$ (14)

and the approximation of $y(t_{n+1})$ is given by $y_{n+1}=y_{n,M}$, since
$c_{M}=1$. A stable way of computing the matrix coefficients $a_{ij}(\Delta
t\Lambda(y_{n}))$ is proposed in Appendix A.

Remark that the first equation of 14, $i=0$, is simply $y_{n,0}=y_{n}$.
Therefore, system 14 can be written in compact form as a $N(M+1)\times N(M+1)$
nonlinear system

$\bm{C}(\bm{y}_{n})\coloneqq\left(I-\Delta tA(\Delta
t\Lambda(y_{n,0}))\bm{G}\right)(\bm{y}_{n})=\mathbf{1}\otimes y_{n},$ (15)

where $\otimes$ is the Kronecker product, $\mathbf{1}\in\mathbb{R}^{M+1}$ a
vector of ones, and $I\in\mathbb{R}^{N(M+1)\times N(M+1)}$ is the identity
matrix. The block matrix $A(\Delta
t\Lambda(y_{n,0}))\in\mathbb{R}^{N(M+1)\times N(M+1)}$, with
$(M+1)\times(M+1)$ blocks, and the nonlinear term
$\bm{G}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{N(M+1)}\rightarrow\mathbb{R}^{N(M+1)}$, are given by

$A(z)=\begin{pmatrix}0&0&\cdots&0\\\ 0&a_{11}(z)&\cdots&a_{1M}(z)\\\
\vdots&\vdots&\ddots&\vdots\\\
0&a_{M1}(z)&\cdots&a_{MM}(z)\end{pmatrix},\qquad\bm{G}(\bm{y}_{n})=\begin{pmatrix}g(y_{n,0})\\\
g(y_{n,1})\\\ \vdots\\\ g(y_{n,M}))\end{pmatrix},$ (16)

where in the definition 8 of $g$ we also replace $y_{n}$ with $y_{n,0}$.
Indeed, the main advantage of including $y_{n,0}$ into $\bm{y}_{n}$ is that
due to the first line of 15 we can interchange $y_{n,0}$ and $y_{n}$. In this
way the definition of $\bm{C}$ doesn’t directly depend on $y_{n}$ and can be
used for all time steps.

A few comments are warranted regarding the interconnections between 15,
collocation methods and exponential Runge–Kutta methods. First, for $z=0$, 13
reduce to the weights of a standard collocation method [21, Section II.7].
Here $z=\Delta t\Lambda(y_{n,0})$ is a nonzero block diagonal matrix 6b;
however, the first two diagonal blocks of $\Lambda(y_{n,0})$ are indeed zero
matrices. Thus, for the first two components $V_{m},\bm{w}_{a}$ of $y$, 15
simplifies a collocation scheme, which motivates our choice of Radau IIA nodes
for $c_{1},\ldots,c_{M}$. In contrast, the last diagonal block
$\Lambda_{g}(V_{m})$ of $\Lambda(y_{n,0})$ is nonzero, thus for the last
component $\bm{w}_{g}$ of $y$, 15 is an exponential Runge–Kutta method of
collocation type [23]. Recall that $\Lambda_{g}(V_{m})$ is diagonal, therefore
its exponentials are inexpensive to compute. We conclude that, for the
monodomain equation, method 15 is an hybrid between collocation schemes and
exponential Runge–Kutta methods.

####  3.1.3 Preconditioned Picard iteration

The monolithic $N(M+1)\times N(M+1)$ nonlinear system 15 is expensive to
solve. Furthermore, for some initial guesses, the Newton method might struggle
to converge. Spectral deferred corrections are alternative iterative methods
for solving 15. These methods were introduced in [14], where they have been
derived by defining a correction equation, which is iteratively solved and
improves the current solution (see also [6, 5, 29, 20, 43, 38, 44, 7]). More
recently, the same methods were derived by solving a preconditioned Picard
iteration [3, 4, 35, 26, 34, 33]. For the purpose of this paper, we follow the
latter approach, due to its convenient notation for the parallel-in-time
extension.

The preconditioned Picard iteration for 15 is given by

$\bm{P}(\bm{y}_{n}^{k+1})=(\bm{P}-\bm{C})(\bm{y}_{n}^{k})+\mathbf{1}\otimes
y_{n},$ (17)

where $\bm{P}$ is the nonlinear preconditioner operator and as initial guess
we choose $\bm{y}_{n}^{0}=\mathbf{1}\otimes y_{n}$. In the remainder of this
section we define the preconditioner $\bm{P}$ for the monodomain equation.

From 10, we have

$\displaystyle y(t_{n}+c_{i}\Delta t)$ $\displaystyle=y_{n}+\Delta
t\int_{0}^{c_{i}}e^{(c_{i}-s)\Delta t\Lambda(y_{n})}g(y(t_{n}+s\Delta
t))\operatorname{d\\!}s$ (18) $\displaystyle=y_{n}+\Delta
t\sum_{j=1}^{i}e^{(c_{i}-c_{j})\Delta
t\Lambda(y_{n})}\int_{c_{j-1}}^{c_{j}}e^{(c_{j}-s)\Delta
t\Lambda(y_{n})}g(y(t_{n}+s\Delta t))\operatorname{d\\!}s.$ (19)

Let $d_{i}=c_{i}-c_{i-1}$ for $i=1,\ldots,M$. Since $e^{(c_{i}-c_{j})\Delta
t\Lambda(y_{n})}=1+\mathcal{O}(\Delta t)$, these terms can be dropped from the
forthcoming computations without affecting accuracy. Therefore, approximations
$\tilde{y}_{n,i}$ to $y(t_{n}+c_{i}\Delta t)$, for $i=0,\ldots,M$, are
obtained as

$\displaystyle\tilde{y}_{n,i}$ $\displaystyle=y_{n}+\Delta
t\sum_{j=1}^{i}\int_{c_{j-1}}^{c_{j}}e^{(c_{j}-s)\Delta
t\Lambda(y_{n})}(f_{I}(\tilde{y}_{n,j})+f_{E}(\tilde{y}_{n,j-1})+f_{e}(\tilde{y}_{n,j-1})+\Lambda(y_{n})(y_{n}-\tilde{y}_{n,j-1}))\operatorname{d\\!}s$
(20) $\displaystyle=y_{n}+\Delta t\sum_{j=1}^{i}d_{j}\varphi_{1}(d_{j}\Delta
t\Lambda(y_{n}))(f_{I}(\tilde{y}_{n,j})+f_{E}(\tilde{y}_{n,j-1})+f_{e}(\tilde{y}_{n,j-1})+\Lambda(y_{n})(y_{n}-\tilde{y}_{n,j-1})),$

with

$\varphi_{1}(z)=\int_{0}^{1}e^{(1-r)z}\operatorname{d\\!}r=\frac{e^{z}-1}{z}.$
(21)

The first two diagonal blocks of $\Lambda(y_{n})$ are zero and
$\varphi_{1}(0)=1$, therefore, from 6a,

$\varphi_{1}(d_{j}\Delta
t\Lambda(y_{n}))(f_{I}(\tilde{y}_{n,j})+f_{E}(\tilde{y}_{n,j-1}))=f_{I}(\tilde{y}_{n,j})+f_{E}(\tilde{y}_{n,j-1})$
(22)

and 20 is

$\tilde{y}_{n,i}=y_{n}+\Delta
t\sum_{j=1}^{i}d_{j}(f_{I}(\tilde{y}_{n,j})+f_{E}(\tilde{y}_{n,j-1})+\varphi_{1}(d_{j}\Delta
t\Lambda(y_{n}))(f_{e}(\tilde{y}_{n,j-1})+\Lambda(y_{n})(y_{n}-\tilde{y}_{n,j-1}))).$
(23)

As previously mentioned, neglecting terms
$e^{(c_{i}-c_{j})\Lambda(y_{n})\Delta t}=1+\mathcal{O}(\Delta t)$ do not
impact accuracy; indeed, the introduced error is comparable, asymptotically,
to the one of the approximated integrals. The reason for their removal is to
improve computational efficiency, since they would require additional
$\mathcal{O}(M^{2})$ diagonal matrix exponential computations. Furthermore, 23
corresponds to the popular IMEX Rush–Larsen method in cardiac
electrophysiology, where the diffusion is integrated implicitly, the stiff
gating variables exponentially, and the remaining nonstiff terms explicitly
[32, 27].

Overall, 23 is significantly cheaper to compute than 14, indeed it is a
diagonally implicit method where the implicit term is linear. Therefore, from
23, we define the preconditioning operator $\bm{P}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{N(M+1)}\rightarrow\mathbb{R}^{N(M+1)}$ as
$\bm{P}(\bm{y}_{n})=(\bm{P}_{0}(\bm{y}_{n}),\ldots,\bm{P}_{M}(\bm{y}_{n}))$,
with

$\displaystyle\bm{P}_{i}(\bm{y}_{n})=y_{n,i}-\Delta t\sum_{j=1}^{i}d_{j}$
$\displaystyle\left(f_{I}(y_{n,j})+f_{E}(y_{n,j-1})\right.$ (24)
$\displaystyle\left.\quad+\varphi_{1}(d_{j}\Delta
t\Lambda(y_{n}))(f_{e}(y_{n,j-1})+\Lambda(y_{n})(y_{n}-y_{n,j-1}))\right)$

for $i=0,\ldots,M$. Remark that for the $f_{I},f_{E}$ terms this is the
standard preconditioner used in SI-SDC methods [29]. However, here we also
have exponential integration of the $f_{e}$ term.

With a close look to 17 and 24, we notice that $y_{n,0}^{k}=y_{n}$ for all
$k$. Hence, for computational efficiency, all terms in 24 depending on
$y_{n,0}$ can be neglected, since they appear on both sides of 17.

###  3.2 Parallel-in-time method

From 17, the derivation of a PinT method based on the PFASST methodology [15,
16] is obtained following the standard procedure based on the Full
Approximation Scheme (FAS) [4, 3, 33, 34, 35], which we recall here for
completeness.

A sequence of $P$ consecutive time steps of the ERK problem 15 can be written
as a $N(M+1)P\times N(M+1)P$ nonlinear system

$\displaystyle\begin{pmatrix}\bm{C}&&&\\\ -\bm{H}&\bm{C}&&\\\
&\ddots&\ddots&\\\ &&-\bm{H}&\bm{C}\end{pmatrix}\begin{pmatrix}\bm{y}_{0}\\\
\bm{y}_{1}\\\ \vdots\\\
\bm{y}_{P-1}\end{pmatrix}=\begin{pmatrix}\mathbf{1}\otimes y_{0}\\\ \bm{0}\\\
\vdots\\\ \bm{0}\end{pmatrix},$ (25)

with $\bm{H}=\bm{N}\otimes I_{N}$ the matrix taking the last node value at the
previous step to be used as initial value in the current step, hence
$\bm{N}\in\mathbb{R}^{(M+1)\times(M+1)}$ has ones in the last column and zeros
elsewhere.

Let $\bm{z}=(\bm{y}_{0},\ldots,\bm{y}_{P-1})\in\mathbb{R}^{N(M+1)P}$ and
$\bm{C}_{P}(\bm{z})=(\bm{C}(\bm{y}_{0}),\ldots,\bm{C}(\bm{y}_{P-1}))^{\top}$,
then 25 can be written as

$\bm{D}(\bm{z})\coloneqq(\bm{C}_{P}-\bm{E}\otimes\bm{H})(\bm{z})=\bm{b},$ (26)

with $\bm{E}\in\mathbb{R}^{P\times P}$ having ones in the first lower off-
diagonal and zeros elsewhere, and $\bm{b}=(\mathbf{1}\otimes
y_{0},\bm{0},\ldots,\bm{0})\in\mathbb{R}^{N(M+1)P}$. In order to solve 26 with
a preconditioned Picard iteration, hence

$\bm{Q}(\bm{z}^{k+1})=(\bm{Q}-\bm{D})(\bm{z}^{k})+\bm{b},$ (27)

there are two types of preconditioners $\bm{Q}$222Throughout this section, the
global ($P$-steps) version of local ($1$-step) operators is denoted by the
subsequent letter; for instance, the global version of $\bm{C}$ is denoted
$\bm{D}$, and the one of $\bm{P}$ is $\bm{Q}$.:

$\displaystyle\bm{Q}^{seq}(\bm{z})$
$\displaystyle=(\bm{P}_{P}-\bm{E}\otimes\bm{H})(\bm{z}),$ (28)
$\displaystyle\bm{Q}^{par}(\bm{z})$ $\displaystyle=\bm{P}_{P}(\bm{z}),$

with
$\bm{P}_{P}(\bm{z})=(\bm{P}(\bm{y}_{0}),\ldots,\bm{P}(\bm{y}_{P-1}))^{\top}$.
Notice as $\bm{Q}^{par}$ can be applied in parallel across the time steps,
while $\bm{Q}^{seq}$ propagates sequentially the solution form one step to the
next.

The PFASST parallel-in-time method employs the full approximation scheme (FAS)
[42] to solve 27 on a multilevel hierarchy, employing the sequential
preconditioner $\bm{Q}^{seq}$ at the coarsest level and the parallel one
$\bm{Q}^{par}$ at any other level. When going from one level to the next,
coarsening can happen in both space and time. Space coarsening is usually done
by decreasing the mesh resolution and/or the polynomial order (as in
$hp$-refinement). Time coarsening is done by keeping the step size $\Delta t$
constant but decreasing the number $M$ of collocation nodes. For the purpose
of this paper, we focus on time coarsening only.

Let $M_{L}<M_{L-1}<\ldots<M_{1}$ be a sequence of collocation nodes’ numbers
and $\\{c_{i}^{l}\\}_{i=1}^{M_{l}}$ the corresponding Radau IIA collocation
nodes, for $l=1,\ldots,L$. Let $\bm{D}_{l}$, and
$\bm{Q}_{l}^{seq},\bm{Q}_{l}^{par}$ be as in 26 and 28, respectively, but for
$M=M_{l}$ collocation nodes. Similarly, let
$\bm{y}_{n}^{l}=(y_{n,0}^{l},\ldots,y_{n,M_{l}}^{l})\in\mathbb{R}^{N(M_{l}+1)}$
and
$\bm{z}^{l}=(\bm{y}_{0}^{l},\ldots,\bm{y}_{P-1}^{l})\in\mathbb{R}^{N(M_{l}+1)P}$.

Let $\bm{I}_{l}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{N(M_{l}+1)}\times\mathbb{R}\rightarrow\mathbb{R}^{N}$ be
the Lagrange interpolator at nodes $c_{1}^{l},\ldots,c_{M_{l}}^{l}$ defined by

$\bm{I}_{l}(\bm{y}_{n}^{l},s)=\sum_{j=1}^{M_{l}}\ell_{j}^{l}(s)y_{n,j}^{l},\qquad\ell_{j}^{l}(s)=\prod_{\begin{subarray}{c}k=1\\\
k\neq j\end{subarray}}^{M_{l}}\frac{s-c_{k}^{l}}{c_{j}^{l}-c_{k}^{l}}.$ (29)

The restriction operator $\bm{R}_{l}^{l+1}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{N(M_{l}+1)}\rightarrow\mathbb{R}^{N(M_{l+1}+1)}$, from
level $l$ to level $l+1$, is linear and we represent it by a matrix,
implicitly defined by

$\bm{R}_{l}^{l+1}\bm{y}_{n}^{l}\coloneqq(y_{n,0}^{l},\bm{I}_{l}(\bm{y}_{n}^{l},c_{1}^{l+1}),\ldots,\bm{I}_{l}(\bm{y}_{n}^{l},c_{M_{l+1}}^{l+1}))^{\top}.$
(30)

The prolongation operator $\bm{T}_{l+1}^{l}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{N(M_{l+1}+1)}\rightarrow\mathbb{R}^{N(M_{l}+1)}$ is
defined as

$\bm{T}_{l+1}^{l}\bm{y}_{n}^{l+1}\coloneqq(y_{n,0}^{l+1},\bm{I}_{l+1}(\bm{y}_{n}^{l+1},c_{1}^{l}),\ldots,\bm{I}_{l+1}(\bm{y}_{n}^{l+1},c_{M_{l}}^{l}))^{\top}.$
(31)

The global restriction operator $\bm{S}_{l}^{l+1}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{N(M_{l}+1)P}\rightarrow\mathbb{R}^{N(M_{l+1}+1)P}$ is
simply $\bm{S}_{l}^{l+1}=I_{P}\otimes\bm{R}_{l}^{l+1}$. Similarly, the global
prolongation operator is $\bm{U}_{l+1}^{l}=I_{P}\otimes\bm{T}_{l+1}^{l}$.

In Algorithm 1 we provide a pseudo code for the multilevel and parallel HSDC
method. Therein, the Burn-in method is an initialization procedure described
in [15]. We stress that Algorithm 1 is a rough description of the true
computational flow of HSDC, which employs the implementation available in the
pySDC [35] library and is slightly different. For instance, in the true code,
if all the residuals up to some step $n$ are below the tolerance, then
iterations at those steps stop and only the remaining steps continue
iterating.

1:function HSDC($y_{0}$, $\Delta t$, $P$, $\\{M_{l}\\}_{l=1}^{L}$, $K$)

2: $\bm{z}_{1}^{0}=$ Burn-in($y_{0}$, $\Delta t$, $P$,
$\\{M_{l}\\}_{l=1}^{L}$)

3: $\bm{b}_{1}=(\mathbf{1}\otimes y_{0},\bm{0},\ldots,\bm{0})$ and
$\bm{b}_{l}=\bm{S}_{l-1}^{l}\bm{b}_{l-1}$, $l=2,\ldots,L$.

4: $k=0$, converged = False

5: while not converged and $k<K$ do

6: $\bm{z}_{1}^{k+1/2}=\bm{z}_{1}^{k}$, $\bm{\tau}_{1}^{k}=\bm{0}$

7: for $l=2,\ldots,L-1$ do

8:
$\bm{\tau}_{l}^{k}=\bm{D}_{l}(\bm{S}_{l-1}^{l}\bm{z}_{l-1}^{k+1/2})-\bm{S}_{l-1}^{l}\bm{D}_{l-1}(\bm{z}_{l-1}^{k+1/2})+\bm{S}_{l-1}^{l}\bm{\tau}_{l-1}^{k}$
$\triangleright$ Residual correction

9: Solve
$\bm{Q}_{l}^{par}(\bm{z}_{l}^{k+1/2})=(\bm{Q}_{l}^{par}-\bm{D}_{l})(\bm{S}_{l-1}^{l}\bm{z}_{l-1}^{k+1/2})+\bm{b}_{l}+\bm{\tau}_{l}^{k}$

10: end for

11:
$\bm{\tau}_{L}^{k}=\bm{D}_{L}(\bm{S}_{L-1}^{L}\bm{z}_{L-1}^{k+1/2})-\bm{S}_{L-1}^{L}\bm{D}_{L-1}(\bm{z}_{L-1}^{k+1/2})+\bm{S}_{L-1}^{L}\bm{\tau}_{L-1}^{k}$

12: Solve
$\bm{Q}_{L}^{seq}(\bm{z}_{L}^{k+1})=(\bm{Q}_{L}^{seq}-\bm{D}_{L})(\bm{S}_{L-1}^{L}\bm{z}_{L-1}^{k+1/2})+\bm{b}_{L}+\bm{\tau}_{L}^{k}$

13: for $l=L-1,\ldots,1$ do

14:
$\bm{c}_{l}^{k+1/2}=\bm{T}_{l+1}^{l}(\bm{z}_{l+1}^{k+1}-\bm{S}_{l}^{l+1}\bm{z}_{l}^{k+1/2})$
$\triangleright$ Coarse grid correction

15: Solve
$\bm{Q}_{l}^{par}(\bm{z}_{l}^{k+1})=(\bm{Q}_{l}^{par}-\bm{D}_{l})(\bm{z}_{l}^{k+1/2}+\bm{c}_{l}^{k+1/2})+\bm{b}_{l}+\bm{\tau}_{l}^{k}$

16: end for

17: $k\leftarrow k+1$

18:
$(\bm{r}_{0}^{k},\ldots,\bm{r}_{P-1}^{k})\coloneqq\bm{D}_{1}(\bm{z}^{k}_{1})-\bm{b}_{1}$
$\triangleright$ Residuals

19: $(\bm{y}_{0}^{k},\ldots,\bm{y}_{P-1}^{k})=\bm{z}^{k}_{1}$ $\triangleright$
Solutions

20: converged = True if $\|\bm{r}_{n}^{k}\|<tol\cdot\|\bm{y}_{n}^{k}\|$ for
$n=0,\ldots,P-1$, else False.

21: end while

22: $y_{P}=y_{P-1,M}$

23: return $y_{P}$

24:end function

Algorithm 1 $P$ steps of the monodomain equation in parallel-in-time.

Remark that the HSDC method is fully defined by the number of levels $L$, the
numbers of collocation nodes $M_{1},\ldots,M_{L}$ at each level, and the
number of processors $P$. To alleviate the notation, henceforth the number of
collocation nodes is noted $M=M_{1},\ldots,M_{L}$. For instance, a two-level
HSDC method with $M_{1}=6$ and $M_{2}=3$ is simply denoted by $M=6,3$.

##  4 Stability analysis

In this section, we begin by numerically assessing the absolute stability of
the HSDC method. Subsequently, we demonstrate that the plain SDC method
exhibits instability when applied to the monodomain equation, thereby
confirming the necessity for HSDC.

###  4.1 Stability of HSDC

We define the linear scalar test equation

$y^{\prime}=\lambda_{I}y+\lambda_{E}y+\lambda_{e}y,\qquad y(0)=1,$ (32)

where $\lambda_{I},\lambda_{E},\lambda_{e}$ represent $f_{I},f_{E},f_{e}$,
respectively. We will verify stability of HSDC after one and also after
multiple time steps, due to the parallel-in-time approach. Therefore, we
define the stability function as

$R_{P}(\lambda_{I},\lambda_{E},\lambda_{e})=y_{P},$ (33)

where $y_{P}$ is the approximation of 32 given by HSDC (Algorithm 1) after $P$
steps of size $\Delta t=1$. The number of levels $L$ and the collocation nodes
$M$ are left out of notation and specified later. Moreover, if $P>1$ we
suppose that steps are computed in parallel, hence with one call to Algorithm
1. The stability domain of the method is defined as

$\mathcal{S}_{P}=\\{\lambda_{I},\lambda_{e},\lambda_{E}\in\mathbb{R}\,\mathrel{\mathop{\mathchar
58\relax}}\,|R_{P}(\lambda_{I},\lambda_{E},\lambda_{e})|\leq 1\\}.$ (34)

Obviously, the exact solution satisfies $|y(P)|\leq 1$ for all
$\lambda_{I},\lambda_{E},\lambda_{e}\leq 0$. For HSDC, we cannot expect
$|y_{P}|\leq 1$ for all $\lambda_{E}\leq 0$ due to the explicit step on the
$f_{E}$ term. However, we don’t expect any lower bound on
$\lambda_{I},\lambda_{e}\leq 0$, since $f_{I},f_{e}$ are integrated implicitly
and exponentially, respectively.

To verify stability of HSDC, we plot
$|R_{P}(\lambda_{I},\lambda_{E},\lambda_{e})|$ for fixed $\lambda_{E}\leq 0$
and varying $\lambda_{I},\lambda_{e}\leq 0$. First, we verify stability in a
nonstiff regime, hence we set $\lambda_{E}=0$ and let
$(\lambda_{I},\lambda_{e})\in[-2,0]\times[-2,0]$. For a single-level HSDC
($L=1$), with $M=6$ collocation nodes, one time step ($P=1$), fixed relative
tolerance $tol=10^{-10}$ and different values of maximum iterations $K=1,2,3$
(cf. Algorithm 1) we display the results in Figs. 1(a), 1(b) and 1(c). Results
for a two-level HSDC with $M=6,3$ are displayed in Figs. 1(d), 1(e) and 1(f).
We consider as well a two-level HSDC with $P=4$ parallel steps, results are
shown in Figs. 1(g), 1(h) and 1(i). We see as in all cases
$|R_{P}(\lambda_{I},\lambda_{E},\lambda_{e})|\leq 1$ and the methods are
therefore stable. Furthermore, when $P=1$ (Figs. 1(a), 1(b), 1(c), 1(d), 1(e)
and 1(f)) the method converges in very few iterations to the same stability
domain, indeed for such small values of $\lambda_{I},\lambda_{E},\lambda_{e}$
we are in a convergence regime. For $P=4$ parallel steps (Figs. 1(g), 1(h) and
1(i)) the scheme converges in few iterations only for those
$\lambda_{I},\lambda_{e}$ in a neighborhood of $(0,0)$, but takes longer to
converge otherwise. Nevertheless, it is stable for all checked values of
$\lambda_{I},\lambda_{e}$.

(a) $L=1$, $P=1$, $K=1$.

(b) $L=1$, $P=1$, $K=2$.

(c) $L=1$, $P=1$, $K=3$.

(d) $L=2$, $P=1$, $K=1$.

(e) $L=2$, $P=1$, $K=2$.

(f) $L=2$, $P=1$, $K=3$.

(g) $L=2$, $P=4$, $K=1$.

(h) $L=2$, $P=4$, $K=2$.

(i) $L=2$, $P=4$, $K=3$.

Figure 1: Stability domain for HSDC in a nonstiff regime, with
$\lambda_{E}=0$, $(\lambda_{e},\lambda_{I})\in[-2,0]\times[-2,0]$ and varying
levels $L$, parallel steps $P$, and maximum iterations $K$. The number of
collocation nodes is $M=6$ for $L=1$ and $M=6,3$ for $L=2$.

In Fig. 2 we perform similar experiments but in a stiff regime. Hence we fix
$\lambda_{E}=-2$ and let
$(\lambda_{I},\lambda_{e})\in[-1000,0]\times[-1000,0]$. We consider the same
combinations of levels $L$, collocation nodes $M$, and processors $P$.
However, for the maximum number of iterations we consider $K=1,2,100$, hence
in the last case we let the algorithm iterate until convergence. We observe
that at convergence ($K=100$) the serial methods ($P=1$) have a very similar
stability domain, compare Figs. 2(c) and 2(f). Interestingly, the stability
domain of the parallel method at convergence, hence Fig. 2(i), has a pattern
similar to the serial methods but with much smaller absolute values. Finally,
note that the HSDC scheme is again stable in all considered cases.

(a) $L=1$, $P=1$, $K=1$.

(b) $L=1$, $P=1$, $K=2$.

(c) $L=1$, $P=1$, $K=100$.

(d) $L=2$, $P=1$, $K=1$.

(e) $L=2$, $P=1$, $K=2$.

(f) $L=2$, $P=1$, $K=100$.

(g) $L=2$, $P=4$, $K=1$.

(h) $L=2$, $P=4$, $K=2$.

(i) $L=2$, $P=4$, $K=100$.

Figure 2: Stability domain for HSDC in a stiff regime, with $\lambda_{E}=-2$,
$(\lambda_{e},\lambda_{I})\in[-1000,0]\times[-1000,0]$ and varying levels $L$,
parallel steps $P$, and maximum iterations $K$. The number of collocation
nodes is $M=6$ for $L=1$ and $M=6,3$ for $L=2$.

###  4.2 Instabilities of SDC for the monodomain equation

We conclude this section by displaying the stability domain of a naive PinT
method for the monodomain equation. Such method would be obtained by using the
standard SDC method in combination with a preconditioner based on the popular
Rush–Larsen method for the monodomain equation [32], hence the one defined in
24. Moreover, recall that the plain SDC method is obtained by replacing
$\Lambda(y_{n})$ in Section 3.1.2 with $0$.

We show the stability domain of this SDC scheme in Fig. 3, where we see that
already after one iteration the method becomes unstable for
$\lambda_{e}<\lambda_{I}$. Therefore, the employment of exponential
Runge–Kutta methods, as in HSDC, is crucial for stability.

(a) $L=1$, $P=1$, $K=1$.

(b) $L=2$, $P=1$, $K=1$.

(c) $L=2$, $P=4$, $K=1$.

Figure 3: Stability domain of a naive SDC method for the monodomain equation,
with $\lambda_{E}=-2$, $(\lambda_{e},\lambda_{I})\in[-1000,0]\times[-1000,0]$
and varying levels $L$ and parallel steps $P$. We fix $K=1$ and the number of
collocation nodes is $M=6$ for $L=1$ and $M=6,3$ for $L=2$. White regions
represent values above $1$, hence instability.

##  5 Numerical experiments

In this section the accuracy, stability and robustness properties of the HSDC
method are assessed numerically. In particular, we first verify the
convergence properties of HSDC for single- and multi-level cases and as well
in parallel. Second, we investigate the impact of the ionic model and its
stiffness on the number of iterations taken by HSDC. Third, we investigate
stability and iteration number of HSDC in a realistic setting, considering
different step sizes, levels and processors. Finally, we display in details
the evolution of the residuals and iterations over time.

##### Model setup.

Throughout the next experiments we will solve the monodomain equation with an
isotropic conductivity
$\sigma=\sigma_{i}\sigma_{e}/(\sigma_{i}+\sigma_{e})\,$\mathrm{mS}\text{\,}{\mathrm{mm}}^{-1}$$,
where $\sigma_{i}=$0.17\text{\,}\mathrm{mS}\text{\,}{\mathrm{mm}}^{-1}$$ and
$\sigma_{e}=$0.62\text{\,}\mathrm{mS}\text{\,}{\mathrm{mm}}^{-1}$$ are the
intracellular and extracellular electric conductivities. The cell membrane
electrical capacitance is set to
$C_{m}=$0.01\text{\,}\mathrm{\SIUnitSymbolMicro
F}\text{\,}{\mathrm{mm}}^{-2}$$ and the surface area-to-volume ratio is
$\chi=$140\text{\,}{\mathrm{mm}}^{-1}$$. These are the same values as in [31].
The physical domain will be
$\Omega=[$0\text{\,}\mathrm{mm}$,$100\text{\,}\mathrm{mm}$]^{d}$. The
dimension $d=1,2$ and the end time $T$ are specified in the experiments below.

In most of the forthcoming simulations, we employ the ten-Tusscher–Panfilov
(TTP) ionic model [41]. It is a very stiff model with $m_{1}=6$ auxiliary
variables and $m_{2}=12$ gating variables. In Section 5.2, where we compare
the effect of the ionic models on the convergence of HSDC, we also consider
the mildly stiff Courtemance–Ramirez–Nattel (CRN) model [12], with $m_{1}=8$
and $m_{2}=12$, and the nonstiff Hodgkin–Huxley (HH) model [25], with
$m_{1}=0$ and $m_{2}=3$. In a typical simulation, with solution away from
equilibrium, the spectral radius of the Jacobian of the ionic models’ right-
hand side is $\rho_{\text{TTP}}\approx 1000$, $\rho_{\text{CRN}}\approx 130$,
and $\rho_{\text{HH}}\approx 40$, for TTP, CRN, and HH, respectively. Note
that the TTP and CRN are realistic models used in practical applications,
while the HH model is mostly used in an academic context. Typical solutions
obtained with the three ionic models are displayed in Fig. 4

For testing the PinT HSDC method we seek for solutions having a fairly
constant, in time, “difficulty”; measured in degrees of stiffness and
nonlinearity, for instance. In particular, we want to avoid convergence to an
equilibrium, in this way we ensure that all processors have to perform
approximately the same amount of work. To do so, we solve monodomain equations
showing a self-sustained wave, as depicted in Fig. 4. There, we also
illustrate the initial value, depending on the ionic model, employed in all
simulations coming next. Since the initial value already contains a self-
sustained wave, there is no need to stimulate the tissue and in all
simulations we set $I_{\text{stim}}(t,\bm{x})\equiv 0$.

The domain $\Omega$ is discretized with uniform mesh sizes of common size (cf.
[9, 31]) $\Delta x_{f}=100/512\,$\mathrm{mm}$\approx$0.2\text{\,}\mathrm{mm}$$
or $\Delta x_{c}=100/256\,$\mathrm{mm}$\approx$0.4\text{\,}\mathrm{mm}$$ and
the Laplacian is approximated with a fourth-order finite difference method.
The time step size $\Delta t$ will usually span a large range of values
specified in each experiment.

(a) TTP, $t=$0\text{\,}\mathrm{ms}$$.

(b) CRN, $t=$0\text{\,}\mathrm{ms}$$.

(c) HH, $t=$0\text{\,}\mathrm{ms}$$.

(d) TTP, $t=$100\text{\,}\mathrm{ms}$$.

(e) CRN, $t=$100\text{\,}\mathrm{ms}$$.

(f) HH, $t=$100\text{\,}\mathrm{ms}$$.

Figure 4: Solutions for different ionic models.

##### Computational setup.

The HSDC method has been implemented in the pySDC [35] library. This library
provides an easy to use and flexible framework for the implementation of
prototypical codes employing SDC and its derivatives, such as MLSDC, PFASST,
and now HSDC. For the particular case of HSDC for the monodomain equation, we
implemented the HSDC sweeper, the tailored preconditioner $\mathbf{P}$ and the
spatial discretization. The linear systems, involving the Laplacian, are
solved via discrete cosine transforms (DCT-II). The ionic models, due to their
high complexity and computational load, are implemented in C++. The code is
found in the official pySDC repository [36]. Numerical experiments are
performed on the Alps(Eiger) cluster at the Swiss National Supercomputing
Centre (CSCS).

###  5.1 Convergence order

In this first experiment we assess the order of convergence of the HSDC
scheme, considering different levels, numbers of collocation nodes, and
processors.

Recall that a collocation method with $M_{1}$ Radau IIA collocation nodes
converges with order of accuracy $p=2M_{1}-1$ [22]. It is also known that the
SDC method with the same nodes converges with order $p=\min(k,2M_{1}-1)$,
where $k$ is the number of iterations [14]. Instead, an exponential
Runge–Kutta collocation method with the same nodes converges with order
$p=M_{1}+1$ (if $M_{1}>1$) [23]. Similarly, ESDC converges with order
$p=\min(k,M_{1}+1)$. Being HSDC an hybrid method between SDC and ESDC, we
expect that it converges with order $p=\min(k,M_{1}+1)$ (we suppose $M_{1}>1$,
hence $M_{1}+1\leq 2M_{1}-1$). However, if the error in the unknowns
integrated with SDC is dominant, then orders of accuracy higher than
$p=\min(k,M_{1}+1)$ can as well be observed.

In general, time integration methods for the monodomain equation are based on
splitting strategies, for efficiency and stability reasons; therefore, those
methods have low order of accuracy due to the splitting error [27]. Higher
order splitting methods are very scarce in the literature, limited to order
$p=4$, and require backward in time steps [8]. For the membrane equation,
hence the monodomain equation without conductivity (3 with $\sigma=0$), higher
order methods are slightly more popular [10, 11, 13] but also limited to order
$p=4$. In contrast, the HSDC method has theoretically arbitrary order of
accuracy. For instance, in this experiment we show that with $M_{1}=6$ nodes
it attains order $p=7$, hence in agreement with the previous considerations.
In some cases it converges with order higher than expected, such as $p=9$ with
$M_{1}=6$. Up to the authors knowledge, the maximum convergence order reached
so far for the monodomain equation was $p=4$.

####  5.1.1 The serial case

Here we perform convergence experiments for the multi-level but serial scheme,
hence $L\geq 1$ and $P=1$. To do so, we solve the monodomain equation 3 in a
two dimensional domain $\Omega=[0,100]\times[0,100]${\mathrm{mm}}^{2}$$, with
$T=$1\text{\,}\mathrm{ms}$$ and the TTP ionic model. We spatially discretize
$\Omega$ with the uniform mesh size $\Delta
x_{c}\approx$0.4\text{\,}\mathrm{mm}$$.

To perform the time convergence experiments, we fix the number of levels $L$,
of collocation nodes $M$, of maximum iterations $K$, and solve the semi-
discrete monodomain equation 5 for different step sizes $\Delta t$. In Fig.
5(a) we display the relative errors against a reference solution with respect
to $\Delta t$, for a single-level HSDC and $M=4$ collocation nodes,
considering different maximum iterations $K$. We see as the order of
convergence increases by one with each additional iteration, attaining the
theoretically maximal order $p=5$. A similar experiment is performed in Fig.
5(b), but with $M=6$ collocation nodes and attaining convergence order $p=7$.
In Fig. 5(c) we show the results for the two-level case, where the order of
convergence increases by two at each iteration. Finally, in Fig. 5(d) the
results for the three-level case are displayed, where the order of convergence
increases by two or even three at each iteration. Note that for the multi-
level cases, we indeed expect to gain more than one order of accuracy per
iteration [26].

(a) $L=1$, $M=4$.

(b) $L=1$, $M=6$.

(c) $L=2$, $M=4,2$.

(d) $L=3$, $M=6,3,1$.

Figure 5: Section 5.1.1. Order of convergence for different levels $L$,
collocation nodes $M$ and maximum number $K$ of iterations.

####  5.1.2 The parallel case

In this experiment perform convergence experiments for the parallel scheme,
hence here $P\geq 1$.

We consider a one dimensional domain $\Omega=[0,100]\,$\mathrm{mm}$$ and set
the end time to $T=$16\text{\,}\mathrm{ms}$$. The domain $\Omega$ is
discretized with the uniform mesh size $\Delta
x_{f}\approx$0.2\text{\,}\mathrm{mm}$$.

Here we will consider a two-level HSDC with $M=6,3$ collocation nodes.
Furthermore, differently from the previous experiment, we do not limit the
number of iterations and let the algorithm iterate till convergence, with a
relative tolerance $tol=5\cdot 10^{-14}$. We display the convergence results
for different processors $P=1,4,16,64$ in Fig. 6(a). All cases converge with
the expected order of convergence $p=7$ and the accuracy is almost independent
from the number of processors. Note that for large step sizes the order of
accuracy is even higher than expected ($p=9$ instead of $p=7$). As previously
mentioned, orders of accuracies higher than $p=M_{1}+1$ are possible if the
error in the $V_{m},w_{a}$ components of $y$ is dominant; indeed, for these
terms the HSDC method corresponds to the SDC method, which is super-convergent
and has order $p=2M_{1}-1$.

In Fig. 6(b) we perform a different convergence experiment where the number of
processors depends on the step size. More specifically, we set $P=T/\Delta t$,
hence all simulation steps are computed simultaneously, disregarding the step
size. Since we consider $\Delta t=2^{-k}$, $k=1,\ldots,7$, then the number of
processors ranges from $P=32$ to $P=2048$. The same orders of accuracies as in
Fig. 6(a) are attained.

(a) Convergence for different $P$.

(b) Convergence for $P$ dependent on $\Delta t$.

Figure 6: Section 5.1.2. Order of convergence in parallel for $L=2$ and
$M=6,3$ collocation nodes. Here HSDC iterates until convergence.

###  5.2 Number of iterations per ionic model

In this experiment we investigate how the choice of the ionic model impacts
the convergence of HSDC. More specifically, we solve the monodomain model on a
two dimensional domain $\Omega=[0,100]\times[0,100]${\mathrm{mm}}^{2}$$ with a
uniform mesh size $\Delta x_{f}\approx$0.2\text{\,}\mathrm{mm}$$ and end time
$T=$1\text{\,}\mathrm{ms}$$. Then, for different HSDC settings and ionic
models, we compare the average number of iterations (over all time steps)
required for convergence, with relative tolerance $tol=5\cdot 10^{-8}$. In
this experiment, we fix $P=1$.

First, in Fig. 7(a) we consider the three ionic models TTP, CRN, and HH, a
single-level HSDC with $M=8$ and vary $\Delta t$. In general, the number of
iterations increases with $\Delta t$, as expected, but remains contained for
the CRN and TTP models. In contrast, for the HH model there is a significant
increase. We believe that this is due to spectral properties: the CRN and TTP
models are mostly dissipative, while HH has oscillatory components [39, 40].

Second, in Fig. 7(b) we fix the step size $\Delta
t=$0.05\text{\,}\mathrm{ms}$$ and let vary the number of collocation nodes $M$
of a single-level HSDC method. Here the average number of iterations increase
as $M$ decreases, since the gap between internal stages increases and thus the
preconditioner $\bm{P}$ becomes less accurate. Note that again CRN and HH are
the less and most demanding models, respectively. Furthermore, for $M>1$ the
TTP model is comparable to CRN.

(a) $M=8$ and varying $\Delta t$.

(b) $\Delta t=$0.05\text{\,}\mathrm{ms}$$ and varying $M$.

Figure 7: Section 5.2. Average number of iterations for different ionic models
and varying $\Delta t$ or collocation nodes $M$. Here a single-level HSDC.

Finally, we consider a two-level HSDC where we fix the number of fine
collocation nodes $M_{1}$ and consider three choices of coarse collocation
nodes $M_{2}$, leading to weak, mild and aggressive coarsening strategies. We
display in Figs. 8(a), 8(b) and 8(c) the number of iterations with respect to
$\Delta t$ and observe as they tend to increase when the number of coarse
collocation nodes $M_{2}$ decreases, as expected. Furthermore, as in the
previous cases, the HH and CRN models are the slowest and fastest to converge,
respectively. For step sizes typically used in practice, hence $\Delta
t\leq$0.1\text{\,}\mathrm{ms}$$, the TTP model also takes few iterations to
converge.

We conclude pointing out that the slow convergence behavior of HSDC for the HH
model in general is not a problem, indeed this model is employed mostly in an
academic context and not in practice. In contrast, the robustness of HSDC for
the CRN and TTP models is much more relevant, since these are very complex
ionic models used in real-life simulations.

(a) $M=8,6$.

(b) $M=8,4$.

(c) $M=8,2$.

Figure 8: Section 5.2. Average number of iterations for different ionic models
and varying $\Delta t$. Here a two-level HSDC with weak, mild and aggressive
coarsening.

###  5.3 Parallel-in-time stability

In Section 4, we showed numerically that for a scalar linear test equation the
HSDC method is stable, for both the single- and multi-level case and serial or
parallel steps. In this experiment we check stability in a more realistic
case, hence solving the monodomain equation 5 on a two dimensional domain
$\Omega=[0,100]\times[0,100]${\mathrm{mm}}^{2}$$ with uniform mesh size
$\Delta x_{f}\approx$0.2\text{\,}\mathrm{mm}$$ and the stiff TTP ionic model.

To carefully asses the stability properties of HSDC, we integrate the problem
employing different settings, by varying $\Delta
t=0.0125,0.025,\ldots,0.2$\mathrm{ms}$$, number of levels $L=2,3,4$ and
parallel steps $P=1,2,\ldots,1024$. For every combination of parameters, we
set the end time to $T=P\Delta t$ and let HSDC iterate till convergence. The
number of iterations required for convergence is an indicator for both
efficiency and stability. In Figs. 9(a), 9(b) and 9(c) we show the results for
$L=2,3,4$, respectively. In each figure, we display the average number of
iterations (over all time steps) with respect to the number of processors $P$
and the shaded areas indicate standard deviation. Moreover, different curves
are associated to a different step size $\Delta t$.

First, we note as the number of iterations increase with the processors $P$,
which is reasonable. For small numbers of processors $P$, the number of
iterations decrease with the additional levels, as expected; indeed, in the
serial case, the multi-level SDC scheme converges in less iterations than
single-level SDC [37]. Nevertheless, this property is not always true in the
parallel case, indeed for large numbers of processors $P$ iterations increase
with the additional levels. The reason comes from the fact that the accuracy
of the Burn-In procedure (cf. Algorithm 1) depends on the coarsest level.
Therefore, if the number of levels is increased by coarsening the coarsest
level, like in this experiment, then the Burn-In procedure becomes less
accurate and provides rougher initial guesses, leading to an increased number
of iterations. This effect is more predominant for long integration intervals,
hence many processors.

In Fig. 9(a) we see that for $P=2$ there is a decrease in the number of
iterations with respect to $P=1$. In general, we expect that iterations
increase with the number of processors. However, for $P=1$ we do not call the
Burn-In procedure, since it is not needed. In contrast, for $P=2$ it is
required and provides initial stage values to both the first and the second
processors. This warmer start in the preconditioned Picard scheme leads to a
decreased number of iterations. The same behavior is not observed in Figs.
9(b) and 9(c) since the Burn-In process is not accurate enough to compensate
for the increased number of iterations required by the additional processor.

Finally, iterations increase also with the step size $\Delta t$. However, a
typical step size for the monodomain equation is $\Delta
t\leq$0.05\text{\,}\mathrm{ms}$$ [9] and here we see that even for $\Delta
t\leq$0.1\text{\,}\mathrm{ms}$$ the number of iterations remains relatively
low compared to the corresponding number of processors.

(a) $L=2$, with $M=8,4$.

(b) $L=3$, with $M=8,4,2$.

(c) $L=4$, with $M=8,4,2,1$.

Figure 9: Section 5.3. Average number of iterations versus number of
processors, for different step sizes and levels.

###  5.4 Number of iterations and residuals

In this final experiment, we show the number of iterations and the residuals
as functions of time. To do so, we discretize again the monodomain equation on
the domain $\Omega=[0,100]\times[0,100]${\mathrm{mm}}^{2}$$ with a mesh size
$\Delta x_{f}\approx$0.2\text{\,}\mathrm{mm}$$ and a time step size $\Delta
t=$0.05\text{\,}\mathrm{ms}$$.

We solve the semi-discrete equation with HSDC in three different settings.
First, a serial and single-level method with $M=8$ collocation nodes. Then, a
serial but two-level HSDC with $M=8,4$ collocation nodes. Finally, a parallel
two-level HSDC with $M=8,4$, and $P=256$. The end time is set to $T=256\Delta
t=$12.8\text{\,}\mathrm{ms}$$ in all cases.

We display in Fig. 10(a) the relative residual
$r_{n}^{k}\coloneqq\|\bm{r}_{n}^{k}\|/\|\bm{y}_{n}^{k}\|$ with respect to time
$t_{n}$ at iteration $k=k_{n}$, where $k_{n}$ is the iterate satisfying the
stopping criteria $\|\bm{r}_{n}^{k_{n}}\|<tol\cdot\|\bm{y}_{n}^{k_{n}}\|$,
i.e. $r_{n}^{k_{n}}<tol$, with $tol=5\cdot 10^{-8}$. In Fig. 10(b) we display
the number of iterations required for convergence, hence $k_{n}$, with respect
to time. Observe as the residual of the single- and two-level serial HSDC
methods oscillates around a constant value, with the single-level method
having smaller residual but performing more iterations. In contrast, the
parallel HSDC method has a residual’s pattern that repeatedly grows in time
and then suddenly drops when an additional iteration is taken. Note as the
parallel method, with $P=256$, resolves the whole time interval $[0,256\Delta
t]$ concurrently, but the number of iterations remains comparable to the
serial methods and increases only slowly in time.

In Figs. 10(c) and 10(d) we display the relative residual $r_{n}^{k}$ with
respect to $k$, for fixed time step $n$ at two different times
$t_{n}=128\Delta t=T/2$ and $t_{n}=256\Delta t=T$. We see as the serial two-
level method is the fastest to converge.

(a) Relative residuals at last iteration $k_{n}$ over time.

(b) Iterations $k_{n}$ needed for convergence over time.

(c) Relative residuals at $t_{n}=$6.4\text{\,}\mathrm{ms}$$.

(d) Relative residuals $t_{n}=$12.8\text{\,}\mathrm{ms}$$.

Figure 10: Section 5.4. Residuals and iterations count over time.

##  6 Conclusion

In this study, we developed a novel variant of the spectral deferred
correction (SDC) method, named Hybrid Spectral Deferred Correction (HSDC).
This method harmoniously blends semi-implicit and exponential SDC techniques,
and is specifically tailored to address the complexities of the monodomain
equation in cardiac electrophysiology. By leveraging the PFASST framework,
HSDC is seamlessly extended to function in a parallel-in-time setting, thereby
enhancing computational efficiency by allowing parallelization across the
temporal direction. This innovative approach not only preserves the essential
stability and robustness needed for such complex simulations but also achieves
an arbitrary order of accuracy.

Our comprehensive numerical experiments underscore the efficacy of HSDC,
particularly highlighting its capability to adeptly manage challenging ionic
models such as the ten-Tusscher–Panfilov. The robustness and accuracy of HSDC
are vital for pushing the boundaries of practical real-time simulations in
cardiac research.

While the primary focus of this paper has been on the development and
validation of the HSDC time discretization scheme, future research should turn
towards code optimization and enhancing spatial discretization. This includes
integrating space-coarsening strategies within the multi-level framework and
exploring the use of the finite element method for spatial discretization.
These advancements will undoubtedly enhance the efficiency and broaden the
applicability of the method.

##  Appendix A Matrix coefficients computation

Here we briefly illustrate how to easily compute the matrix coefficients
$a_{ij}(z)$ from 13. The approach proposed here takes inspiration from [6].

Given any polynomial $p(s)$ of degree $M-1$, the Fornberg algorithm [18] can
be used to write its derivatives $p^{(k)}(0)$ at $s=0$ in terms of its values
at $c_{j}$, $j=1,\ldots,M$:

$p^{(k)}(0)=\sum_{l=1}^{M}w^{k}_{l}p(c_{l}).$ (35)

From a Taylor expansion at $s=0$ and using 35 and
$\ell_{j}(c_{l})=\delta_{jl}$, we have

$\displaystyle\ell_{j}(s)=\sum_{k=1}^{M}\frac{s^{k-1}}{(k-1)!}\ell_{j}^{(k-1)}(0)=\sum_{k=1}^{M}\frac{s^{k-1}}{(k-1)!}\sum_{l=1}^{M}w^{k-1}_{l}\ell_{j}(c_{l})=\sum_{k=1}^{M}\frac{s^{k-1}}{(k-1)!}w^{k-1}_{j}.$
(36)

Plugging 36 into 13 yields

$\displaystyle
a_{ij}(z)=\int_{0}^{c_{i}}e^{(c_{i}-s)z}\ell_{j}(s)\operatorname{d\\!}s=\sum_{k=1}^{M}w_{j}^{k-1}\frac{1}{(k-1)!}\int_{0}^{c_{i}}e^{(c_{i}-s)z}s^{k-1}\operatorname{d\\!}s=\sum_{k=1}^{M}c_{i}^{k}w_{j}^{k-1}\varphi_{k}(c_{i}z)$
(37)

where we performed the change of variables $s=c_{i}r$ and used

$\varphi_{k}(z)\coloneqq\frac{1}{(k-1)!}\int_{0}^{1}e^{(1-r)z}r^{k-1}\operatorname{d\\!}r\quad\text{for
}k\geq 1.$ (38)

Therefore, the matrix coefficients $a_{ij}(z)$ can be computed using 37. The
coefficients $w_{j}^{k}$ are computed with the Fornberg algorithm [18] and
$\varphi_{k}(z)$ by approximating the integral in 38 with a quadrature rule
(we use Gaussian quadrature). The matrix exponential in the integrand of 38 is
easily computable since in our case the matrices are diagonal. Notice that
computing $a_{ij}(z)$ with formula 37 is more stable than using 13, as
numerical evaluation of the highly oscillatory polynomials $\ell_{j}(t)$ is
avoided.

## References

  * [1] M. Bedez, Z. Belhachmi, O. Haeberlé, R. Greget, S. Moussaoui, J.-M. Bouteiller, and S. Bischoff.  A fully parallel in time and space algorithm for simulating the electrical activity of a neural tissue.  Journal of Neuroscience Methods, 257:17–25, 2016.
  * [2] P. Benedusi, M. L. Minion, and R. Krause.  An experimental comparison of a space-time multigrid method with PFASST for a reaction-diffusion problem.  Computers & Mathematics with Applications, 99:162–170, 2021.
  * [3] M. Bolten, D. Moser, and R. Speck.  A multigrid perspective on the parallel full approximation scheme in space and time.  Numerical Linear Algebra with Applications, 24(6):e2110, 2017.
  * [4] M. Bolten, D. Moser, and R. Speck.  Asymptotic convergence of the parallel full approximation scheme in space and time for linear problems.  Numerical Linear Algebra with Applications, 25(6):e2208, 2018.
  * [5] A. Bourlioux, A. T. Layton, and M. L. Minion.  High-order multi-implicit spectral deferred correction methods for problems of reactive flow.  Journal of Computational Physics, 189(2):651–675, 2003.
  * [6] T. Buvoli.  A class of exponential integrators based on spectral deferred correction.  SIAM Journal on Scientific Computing, 42(1):A1–A27, 2020.
  * [7] M. F. Causley and D. C. Seal.  On the convergence of spectral deferred correction methods.  Communications in Applied Mathematics and Computational Science, 14(1):33–64, 2019.
  * [8] J. Cervi and R. J. Spiteri.  High-Order Operator-Splitting Methods for the Bidomain and Monodomain Models.  In D. Boffi, L. F. Pavarino, G. Rozza, S. Scacchi, and C. Vergara, editors, Mathematical and Numerical Modeling of the Cardiovascular System and Applications, volume 16 of SEMA SIMAI Springer Series, pages 23–40. Springer International Publishing, Cham, 2018.
  * [9] R. H. Clayton, O. Bernus, E. M. Cherry, H. Dierckx, F. H. Fenton, L. Mirabella, A. V. Panfilov, F. B. Sachse, G. Seemann, and H. Zhang.  Models of cardiac tissue electrophysiology: Progress, challenges and open questions.  Progress in Biophysics and Molecular Biology, 104(1-3):22–48, 2011.
  * [10] Y. Coudière, C. Douanla-Lontsi, and C. Pierre.  Exponential Adams–Bashforth integrators for stiff ODEs, application to cardiac electrophysiology.  Mathematics and Computers in Simulation, 153:15–34, 2018.
  * [11] Y. Coudière, C. D. Lontsi, and C. Pierre.  Rush-Larsen time-stepping methods of high order for stiff problems in cardiac electrophysiology.  Electronic Transactions on Numerical Analysis, 52:342–357, 2020.
  * [12] M. Courtemanche, R. J. Ramirez, and S. Nattel.  Ionic mechanisms underlying human atrial action potential properties: Insights from a mathematical model.  American Journal of Physiology-Heart and Circulatory Physiology, 275(1):H301–H321, 1998.
  * [13] C. Douanla Lontsi, Y. Coudière, and C. Pierre.  Efficient high order schemes for stiff ODEs in cardiac electrophysiology.  Revue Africaine de Recherche en Informatique et Mathématiques Appliquées, 28, 2018.
  * [14] A. Dutt, L. Greengard, and V. Rokhlin.  Spectral deferred correction methods for ordinary differential equations.  BIT Numerical Mathematics, 40(2):241–266, 2000.
  * [15] M. Emmett and M. L. Minion.  Toward an efficient parallel in time method for partial differential equations.  Communications in Applied Mathematics and Computational Science, 7(1):105–132, 2012.
  * [16] M. Emmett and M. L. Minion.  Efficient Implementation of a Multi-Level Parallel in Time Algorithm.  In J. Erhel, M. J. Gander, L. Halpern, G. Pichot, T. Sassi, and O. Widlund, editors, Domain Decomposition Methods in Science and Engineering XXI, volume 98 of Lecture Notes in Computational Science and Engineering, pages 359–366. Springer International Publishing, Cham, 2014.
  * [17] R. FitzHugh.  Mathematical models of threshold phenomena in the nerve membrane.  The Bulletin of Mathematical Biophysics, 17:257–278, 1955.
  * [18] B. Fornberg.  Generation of Finite Difference Formulas on Arbitrarily Spaced Grids.  Mathematics of Computation, 51(184):699–706, 1988.
  * [19] M. J. Gander and M. Neumüller.  Analysis of a New Space-Time Parallel Multigrid Algorithm for Parabolic Problems.  SIAM Journal on Scientific Computing, 38(4):A2173–A2208, 2016.
  * [20] T. Hagstrom and R. Zhou.  On the spectral deferred correction of splitting methods for initial value problems.  Communications in Applied Mathematics and Computational Science, 1(1):169–205, 2006.
  * [21] E. Hairer, S. P. Nørsett, and G. Wanner.  Solving Ordinary Differential Equations I, volume 8.  Springer-Verlag, Berlin, second edition, 2008.
  * [22] E. Hairer and G. Wanner.  Solving Ordinary Differential Equations II, volume 14.  Springer-Verlag, Berlin, second edition, 2002.
  * [23] M. Hochbruck and A. Ostermann.  Exponential Runge–Kutta methods for parabolic problems.  Applied Numerical Mathematics, 53(2-4):323–339, 2005.
  * [24] M. Hochbruck and A. Ostermann.  Exponential integrators.  Acta Numerica, 19:209–286, 2010.
  * [25] A. L. Hodgkin and A. F. Huxley.  A quantitative description of membrane current and its application to conduction and excitation in nerve.  The Journal of Physiology, 117(4):500–544, 1952.
  * [26] G. Kremling and R. Speck.  Convergence analysis of multi-level spectral deferred corrections.  Communications in Applied Mathematics and Computational Science, 16(2):227–265, 2021.
  * [27] L. P. Lindner, T. Gerach, T. Jahnke, A. Loewe, D. Weiss, and C. Wieners.  Efficient time splitting schemes for the monodomain equation in cardiac electrophysiology.  International Journal for Numerical Methods in Biomedical Engineering, 39(2):e3666, 2023.
  * [28] J. Lions, Y. Maday, and G. Turinici.  Résolution d’EDP par un schéma en temps $\ll$pararéel$\gg$.  Comptes Rendus de l’Académie des Sciences, 332(7):661–668, 2001.
  * [29] M. L. Minion.  Semi-implicit spectral deferred correction methods for ordinary differential equations.  Communications in Mathematical Sciences, 1(3):471–500, 2003.
  * [30] J. Nagumo, S. Arimoto, and S. Yoshizawa.  An Active Pulse Transmission Line Simulating Nerve Axon.  Proceedings of the IRE, 50(10):2061–2070, 1962.
  * [31] S. A. Niederer, E. Kerfoot, A. P. Benson, M. O. Bernabeu, O. Bernus, C. Bradley, E. M. Cherry, R. Clayton, F. H. Fenton, A. Garny, E. Heidenreich, S. Land, M. Maleckar, P. Pathmanathan, G. Plank, J. F. Rodríguez, I. Roy, F. B. Sachse, G. Seemann, O. Skavhaug, and N. P. Smith.  Verification of cardiac tissue electrophysiology simulators using an N-version benchmark.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 369:4331–4351, 2011.
  * [32] S. Rush and H. Larsen.  A Practical Algorithm for Solving Dynamic Membrane Equations.  IEEE Transactions on Biomedical Engineering, (4):389–392, 1978.
  * [33] O. Sander, R. Schöbel, and R. Speck.  The Parallel Full Approximation Scheme in Space and Time for a Parabolic Finite Element Problem.  In Domain Decomposition Methods in Science and Engineering XXVI, pages 531–538. Brenner, Susanne C. and Chung, Eric and Klawonn, Axel and Kwok, Felix and Xu, Jinchao and Zou, Jun, 2022.
  * [34] R. Schöbel and R. Speck.  PFASST-ER: Combining the parallel full approximation scheme in space and time with parallelization across the method.  Computing and Visualization in Science, 23:1–12, 2020.
  * [35] R. Speck.  Algorithm 997: PysDC—prototyping spectral deferred corrections.  ACM Transactions on Mathematical Software, 45(3):1–23, 2019.
  * [36] R. Speck, T. Lunet, T. Baumann, L. Wimmer, and I. Akramov.  Parallel-in-time/pySDC (5.4.3).  Zenodo, 2024.  URL https://doi.org/10.5281/zenodo.10886512.
  * [37] R. Speck, D. Ruprecht, M. Emmett, M. Minion, M. Bolten, and R. Krause.  A multi-level spectral deferred correction method.  BIT Numerical Mathematics, 55(3):843–867, 2015.
  * [38] R. Speck, D. Ruprecht, M. Minion, M. Emmett, and R. Krause.  Inexact Spectral Deferred Corrections.  In T. Dickopf, M. J. Gander, L. Halpern, R. Krause, and L. F. Pavarino, editors, Domain Decomposition Methods in Science and Engineering XXII, volume 104 of Lecture Notes in Computational Science and Engineering, pages 389–396. Springer International Publishing, Cham, 2016.
  * [39] R. J. Spiteri and R. C. Dean.  Stiffness analysis of cardiac electrophysiological models.  Annals of Biomedical Engineering, 38(12):3592–3604, 2010.
  * [40] R. J. Spiteri and R. C. Dean.  Erratum: Stiffness analysis of cardiac electrophysiological models.  Annals of Biomedical Engineering, 40(7):1622–1625, 2012.
  * [41] K. H. W. J. ten Tusscher and A. V. Panfilov.  Alternans and spiral breakup in a human ventricular tissue model.  American Journal of Physiology. Heart and Circulatory Physiology, 291(3):H1088–1100, 2006.
  * [42] U. Trottenberg, C. W. Oosterlee, A. Schüller, A. Brandt, P. Oswald, and K. Stüben.  Multigrid.  Elsevier Academic Press, Amsterdam Heidelberg, transferred to digital print edition, 2007.
  * [43] M. Weiser.  Faster SDC convergence on non-equidistant grids by DIRK sweeps.  BIT Numerical Mathematics, 55(4):1219–1241, 2015.
  * [44] M. Weiser and S. Ghosh.  Theoretically optimal inexact spectral deferred correction methods.  Communications in Applied Mathematics and Computational Science, 13(1):53–86, 2018.


# Yin Yang Convolutional Nets: Image Manifold Extraction by the Analysis of
Opposites

Augusto Seben da Rosa1, Frederico Santos de Oliveira2,
Anderson da Silva Soares3, Arnaldo Candido Junior4

1 Federal University of Technology - Paraná
2 Federal University of Mato Grosso
3 Federal University of Goias
4 São Paulo State University

###### Abstract

Computer vision in general presented several advances such as training
optimizations, new architectures (pure attention, efficient block, vision
language models, generative models, among others). This have improved
performance in several tasks such as classification, and others. However, the
majority of these models focus on modifications that are taking distance from
realistic neuroscientific approaches related to the brain. In this work, we
adopt a more bio-inspired approach and present the Yin Yang Convolutional
Network, an architecture that extracts visual manifold, its blocks are
intended to separate analysis of colors and forms at its initial layers,
simulating occipital lobe’s operations. Our results shows that our
architecture provides State-of-the-Art efficiency among low parameter
architectures in the dataset CIFAR-10. Our first model reached 93.32% test
accuracy, 0.8% more than the older SOTA in this category, while having 150k
less parameters (726k in total). Our second model uses 52k parameters, losing
only 3.86% test accuracy. We also performed an analysis on ImageNet, where we
reached 66.49% validation accuracy with 1.6M parameters. We make the code
publicly available at: https://github.com/NoSavedDATA/YinYang_CNN.

##  1 Introduction

The field of neural computer vision have presented a great advancement, for
instance, improved neural network architectures were proposed. Some of these
architectures include: mobile model families as MobileNet [8], EfficientNet V2
[21] and RegNet [14]; two-branch neural networks for semantic segmentation, as
BiSeNet [25], Deep Dual-Resolution networks [12] and SeaFormer [24]; pure
attention mechanisms applied into image classification, such as ViT [4] and
MaxViT [22]; image generative models like Stable Diffusion [16] and DALL-E-2
[15]; and lastly, vision-language models, as CLIP [13]. The majority of these
architectures focuses on increasing model efficiency by improving the micro-
architecture – there is, by making adjustments relative to the inside of a
network block, as in mobile model families. Some of these architecture also
leverage the potential of CNNs and Transformers into high level computer
vision tasks, as image generation or vision-language models creation.

However, none of the famous modern architectures are trying to reach a more
realistic neuroscientific approach in respect to brain occipital (visual
cortex) mechanisms. In this regard, we base our research on two
neuroscientific findings about function specialized occipital lobe areas, that
is, edge detection that happens on V1 area [5] and color processing at V4 area
[1]. Also, this form of specialization is also observed in the human eye, in
which rod components are related to white and black processing [2], and cone
components to color processing [3].

In this research, we present Yin Yang Convolutional Net (YYNet), a neural
network model which makes adjustments relative to the global scale of the
model. This is performed in the macro-architecture by aggregating blocks (or
single purpose networks) that extracts visual manifolds by doing a separate
analysis of colors and forms from its input. We find that this architecture
provides State-of-the-Art (SOTA) efficiency among low parameter architectures
applied to the low data and image resolution dataset, CIFAR-10, using less
parameters and less training epochs than existing models.

##  2 Related Work

Vision neural network research have grown in quality at a very fast pace. In
this section, we present related architectures and their given tasks. Related
work can be categorized into five approaches: mobile networks; two branch
networks; transformers; vision-language models; generative models.

Some of this approaches are usually divided in three parts: stem, stage and
head. The stem is usually single convolution with stride 2, but may present
optional extra convolutions. The stage contains the main architecture of the
model, that can be divided into blocks with layers, in which each stage block
shares hyperparameters (number of channels and extra hyperparameters) across
all its layers. The head may or not contain convolutions and then it is
followed by average pooling, an optional linear layer and the final
classification linear layer followed by a softmax. This is the case for Mobile
Nets and MaxViT.

###  2.1 Mobile Networks

Mobile network based approaches focus on building efficient models with fewer
parameters. The authors of MobileNet V2 [17], shown in Figure 1, propose
residual inverted bottleneck blocks, a micro-architecture that became more
efficient than the reference model ResNet [6]. The objective of this micro-
architecture is to reduce data dimensionality in a way that the manifold spans
the entire space of lower dimensional sub-spaces. They do so by inserting
inverted bottlenecks at each block, instead of keeping the channels number
constant at repeating blocks as in ResNet. That is, they first expand the
number of channels for a given factor (e.g. 4), then they apply convolutions
in this higher dimension and finally go back into the original dimension,
similar to feed-forward networks in the transformer architecture [23]. On
their architecture, they first increase the channels number with 1x1 kernels,
then use 3x3 depth-wise kernels on the higher dimension and, lastly, they
reduce the dimension back to what it was before with a 1x1 kernel. The reason
to use 1x1 kernels follwed by 3x3 is because the authors of Mobile Net V1 [9]
found it more efficient than directly expanding using 3x3 kernels.

(a) Mobile Net V2

(b) Mobile Net V3

Figure 1: Mobile Net Micro-Architectures [8]

Further, on Mobile Net V3 [8], they improve Mobile Net V2 efficiency by
applying a Squeeze and Excitation [10] mechanism after the 3x3 depth-wise
convolution. We will refer to the block of Mobile Net V3 as the MBConv.

Besides, in EfficientNet V2 [21], the authors have changed the original
architecture of MobileNet V3 into the Fused MBConv and interpolated this new
block with the original MBConv, finding the best parameter configuration given
their applied model search. They have also proposed slight modifications to
the compound scaling method of EfficientNet [20], which consists in adjusting
model depth, channel number and image size to find the best scale inside a
family of models.

###  2.2 Two-branch Neural Networks

Two-branch based models allow the signal to travel along two different paths,
conventionally after a shared stem - except for BiSeNet. This is done to
extract manifold efficiently. An example of different paths for signal
propagation can be seem in BiSeNet [25] and Deep Dual-Resolution Networks
[12], these architectures are built for the purpose of semantic segmentation,
in which one of the branches is shallow and wide, extracting local details,
and the other is deep and narrow, capturing high-level semantics. For
simplicity, we show only the Deep Dual-Resolution Network at Figure 2.

Figure 2: Deep Dual-Resolution Network [12].

Another example of two branch network is SeaFormer [24], that increases the
efficiency by using mobile transformers, from which the embedding of the stem
backbone is processed and fused at multiple steps with mobile transformer
blocks embeddings, as can be seen at Figure 3. Our work is inspired by these
architectures.

Figure 3: SeaFormer [24].

###  2.3 Vision Transformers

Vision Transformers apply transformer encoder blocks to explore the attention
paradigm into vision tasks. On ViT [4], as shown in Figure 4, the authors do
this by dividing the image into multiple embedding patches of equal window
size, in which each patch embedding is cross-attended with attention to all
other patches. They also reserve a patch embedding for classification (usually
referred as “cls” token in language models), which does not come from the
image, and so aggregates global information about other patches. The advantage
of using a pure transformers architecture is that the multiple heads mechanism
turns it possible to apply tensor parallelism [19].

Figure 4: ViT [4].

Then, MaxViT [22] enhances ViT with more efficient attention mechanisms for
the vision paradigm, the authors use Axial Attention for local details and
Grid Attention for global interactions between pixels. They also append a
MBConv at the start of each block.

##  3 Architecture

Yin-Yang Net uses a micro-architcture presented in Figure 5. In each repeating
block, we start with a sub-block of ResNet and then use $n$ sub-blocks of
MBConv from Mobile Net V3. We found this configuration has better accuracy
when training solely with MBConv given our hyperparameters, and it is more
efficient during training than using only ResNet blocks. At the two branch
layers, stride 2 is applied on the last or the first sub-block according to
the micro-architecture type, Yin or Yang, respectively. The single branch
layers use stride 2 on the ResNet block.

Figure 5: Micro-architecture

(a) Micro-arch

(b) Yin

(c) Yang

(d) Single Path

Our work is inspired by two branch architectures. However, our approach
differs from classical two branch networks in 2 aspects. First, instead of
building a shallow and a deep branch for details and semantics extraction,
YYNet uses the same number of layers and channels at blocks on the same level,
but stride 2 at different parts of these layers. Second, in our work, there is
no common stem backbone to the branches, the stem is the focus of this paper,
where different manifolds are analyzed.

The Yin branch has the purpose of form analysis. Yin blocks can use the first
channel of the input or the mean of all channels. We found that both
configurations perform well. For simplicity, when the first channel approach
is in use, the first block receives the red channel in the network input
layer. This way, as there is no other color to extract, it is obliged to the
task of extracting the form manifold. Also, in order to focus on local/higher
scale details, a strategy of later stride 2 is used, meaning that the last
MBConv of each block applies striding.

On the other hand, the Yang branch analyzes colors. With that in mind, as
colors in nearly pixel are generally the same, we use an early stride 2 on
this branch to remove color redundancy and only care about different colors
interactions. That is, the first MBConv on each block applies striding. This
block resembles standard single branch architectures, as its input is the
three RGB channels and early stride 2 is applied.

Then, at the macro-architecture level, as shown in Figure 6, we send the same
input to both these micro-architectures. Their final embeddings have the same
shape, as they have the same number of layers and channels, with the exception
of the first sub-block of each micro-architecture. We then apply a Fusion Gate
mechanism adapted for our embeddings.

Figure 6: Macro-architecture.

We apply a Fusion Gate, similar to SeaFormer and Multimodal Chain-of-Thought
[27]. For this, we use an embedding fusion mechanism at outputs from Yin and
Yang branches, each branch having a different embedding meaning. This is done
to unify both embeddings, or the manifold, as presented in (1).

$SP_{X}=A_{Y}+I_{Y}$ (1)

$X$ represents the input of one network and $Y$ represents an output, $A$
represents the Yang blocks, $I$ represents the Yin blocks, $\odot$ represents
a Hadamard product (elementwise product). We have tested several combinations
of $A$ and $I$ (as described in Section 4). The operations presented in (1)
are those with best performance. On CIFAR-10, this gated fusion yielded a
slightly better accuracy than the concatenation – about less than 0.5%, yet
halving the channels number compared to concatenation.

After this, we send the embedding into Single Path blocks, that consists of a
sub-block of ResNet with stride 2 on the first convolution and then $n$ sub-
blocks of MBConv with no stride 2 (except for CIFAR-10, as described in
Section 4). We use GELU [7] as the activation function of any sub-block. As
the head of our model, we use average pooling, flattening, a linear layer,
GELU, dropout and the final classification linear layer followed by a softmax.

##  4 Experiments and Results

###  4.1 Experiments

We tested multiple fusion approaches to combine Yin and Yang outputs. We
selected the approach presented in Equation 1. However, Table 1 presents other
approaches tested in this work. We performed 3 runs with batch size 512 on
CIFAR-10 for each approach and report the mean and standard deviation.

Table 1: Fusion Approach Formula | Mean | STD
---|---|---
A*(1-I) | 87.61 | 0.09
A*I + A+I | 87.63 | 0.08
A*(1-I) + A-I | 87.81 | 0.19
A*I | 87.87 | 0.23
A*(1-I) + A+I | 87.98 | 0.22
A+I | 88.21 | 0.46

We then reduced the batch size into 64 for the final model on CIFAR-10. We do
this because batch normalization seems to work better in batch sizes on the
range of 50 to 100 [26]. We use a smaller batch size in ImageNet due to
computational constraints.

We used gradient clip of 1 and mixed precision. At 25% of training, we
activate a exponential moving average with a multiplier of 0.1 for the
averaged model parameter and 0.9 for current model parameter. Also, one of the
graphs at [11] shows that there are specific values of weight decay that works
better with specific values of learning rate. We therefore use an adaptive
value for the weight decay. At the end of each epoch, we set the weight decay
to be equal to the learning rate*1.56. The input resolution for CIFAR-10 is
32x32 and ImageNet is 224x224. Further hyperparameters and settings are
provided in Table 2. Specific Hyperparameters adjustments for CIFAR-10 and
ImageNet are presented on Table 3. We designed 3 models for CIFAR-10 and one
for ImageNet.

Table 2: General Hyperparameters and Settings Hyperparameter | YYNet Small | YYNet
---|---|---
Optimizer | AdamW [11] | AdamW
LR Scheduler | One Cycle | One Cycle
Max LR | 1e-2 | 18e-4
Epochs | 40 | 300
Batch Size | 64 | 32
GPU | RTX 2060 | RTX 2080 TI
Dataset | CIFAR-10 | ImageNet
Table 3: Dataset Specific Hyperparameters Hyperparameter | CIFAR-10 | ImageNet
---|---|---
YY Starting Channels | (16, 32, 64) | 16
SP Starting Channels | (16, 32, 64) | 64
Channels added per MBConv | 0 | 2
Extra SP stride 2 | Yes | No
YY Layers | 1 | 1
SP Layers | 1 | 4
YY MBConv per Layer | 3 | 3
SP MBConv per Layer | 2 | 2
Pre-Classification Linear Neurons | 40 | 500

Regarding CIFAR-10, we use a constant channel number over all the sub-blocks,
exploring 3 models variants. They have 1 Yin Yang layer with 3 MBConvs and 1
single branch layer with 2 MBConvs. We apply an extra stride 2 at the first
MBConv of the single branch layer on CIFAR-10 networks. The linear layer
before the classification layer on CIFAR-10 has 40 neurons.

Regarding ImagetNet, we start with 16 channels, then a constant number of 2
channels is added at each MBConv on the Yin and Yang branches. After that, a
fixed number of channels 64 is provided for the first single branch sub-block.
We then continue adding channels after each block. We also use a single layer
and 3 MBConvs for the Yin Yang branches and 2 MBConvs for the single branch
sub-blocks, but in this dataset we use 4 layers of single-branch and no extra
stride 2 is applied. The linear layer before the classification layer has 500
neurons.

###  4.2 Results

We conduct experiments with the small version of YYNet on CIFAR-10, in which
we reach State-of-the-Art (SOTA) at model efficiency for models with few
parameters. These results are provided in Table 4.

Table 4: CIFAR-10 Model | Test Accuracy | Parameters
---|---|---
ExquisiteNetV2 [28] | 92.52 | 890,000
YYNet Small 64 channels (ours) | 93.32 | 726,274
kMobileNet 16ch [18] | 89.81 | 240,000
YYNet Small 32 channels | 91.91 | 191,330
YYNet Small 16 channels | 89.46 | 52,882

We also test a model version on the ImageNet validation dataset. We do not
provide results on the test dataset, since ImageNet team only send results on
the test set when the challenge is open. Comparison with similar size models
are provided in Table 5. Our model uses 6% less parameters than [24] reaching
an validation accuracy only 1.2% smaller. Regarding with MobileNet V3 [8], the
authors do not provide validation accuracy but we assume the validation
accuracy is similar to test accuracy. In this case, our model is considerable
smaller and presents similar efficiency.

Table 5: ImageNet Model | Validation Acc | Test Acc | Parameters
---|---|---|---
MobileNet V3 (Small) [8] | - | 67.4 | 2.5M
SeaFormer (Tiny) [24] | 67.7 | 67.9 | 1.7M
YYNet (ours) | 66.49 | - | 1.6M

##  5 Conclusion and Discussion

In this work, we took inspiration in neuroscience research to model efficient
neural networks. We developed a two branch stem for CNNs intended to analyze
colors and shapes of images separately. Our model reached State-of-the-Art
results on CIFAR-10 considering models with few parameters. We reached 93.32%
test accuracy with 726k parameters, 0.8% more than the older SOTA in this
category, while having close to 150k less parameters. Our model with 52k
parameters, 17 times smaller than ExquisiteNet, lose only 3.86% test accuracy.
We reached 66.49% validation accuracy on ImageNet with 1.6M parameters.

Future work includes parameter search for the Yin Yang network at the ImageNet
dataset. Yin branch would also benefits from this parameter search, as it
applies stride 2 relatively late in this branch, increasing processing cost.

We also plan to investigate if our architecture is useful for other tasks
beyond classification. One example is applying YYNets in generative AIs such
as Stable Diffusion, by changing the latent space currently generated by
U-Nets [16]. Another possible use is combining YYNets shape and color
separation with architectures such as ViT, i.e, adding gray-scale input
patches or queries.

## References

  * [1] M. M. Bannert and A. Bartels.  Human v4 activity patterns predict behavioral performance in imagery of object color.  Journal of Neuroscience, 38(15):3657–3668, 2018.
  * [2] D. Beatty.  Visual pigments and the labile scotopic visual system of fish.  Vision research, 24(11):1563–1573, 1984.
  * [3] B. R. Conway.  Color vision, cones, and color-coding in the cortex.  The neuroscientist, 15(3):274–290, 2009.
  * [4] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale, 2021.
  * [5] J. H. Elder and A. J. Sachs.  Psychophysical receptive fields of edge detection mechanisms.  Vision Research, 44(8):795–813, 2004.
  * [6] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [7] D. Hendrycks and K. Gimpel.  Gaussian error linear units (gelus).  arXiv preprint arXiv:1606.08415, 2016.
  * [8] A. Howard, M. Sandler, G. Chu, L.-C. Chen, B. Chen, M. Tan, W. Wang, Y. Zhu, R. Pang, V. Vasudevan, et al.  Searching for mobilenetv3.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 1314–1324, 2019.
  * [9] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam.  Mobilenets: Efficient convolutional neural networks for mobile vision applications.  arXiv preprint arXiv:1704.04861, 2017.
  * [10] J. Hu, L. Shen, and G. Sun.  Squeeze-and-excitation networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7132–7141, 2018.
  * [11] I. Loshchilov and F. Hutter.  Decoupled weight decay regularization.  In International Conference on Learning Representations, 2019.
  * [12] H. Pan, Y. Hong, W. Sun, and Y. Jia.  Deep dual-resolution networks for real-time and accurate semantic segmentation of traffic scenes.  IEEE Transactions on Intelligent Transportation Systems, 2022.
  * [13] A. Radford, J. W. Kim, C. Hallacy, A. Ramesh, G. Goh, S. Agarwal, G. Sastry, A. Askell, P. Mishkin, J. Clark, et al.  Learning transferable visual models from natural language supervision.  In International conference on machine learning, pages 8748–8763. PMLR, 2021.
  * [14] I. Radosavovic, R. P. Kosaraju, R. Girshick, K. He, and P. Dollar.  Designing network design spaces.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [15] A. Ramesh, P. Dhariwal, A. Nichol, C. Chu, and M. Chen.  Hierarchical text-conditional image generation with clip latents, 2022\.
  * [16] R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer.  High-resolution image synthesis with latent diffusion models.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10684–10695, 2022.
  * [17] M. Sandler, A. Howard, M. Zhu, A. Zhmoginov, and L.-C. Chen.  Mobilenetv2: Inverted residuals and linear bottlenecks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4510–4520, 2018.
  * [18] J. P. Schwarz Schuler, S. Romaní, M. Abdel-nasser, H. Rashwan, and D. Puig.  Grouped pointwise convolutions reduce parameters in convolutional neural networks.  Mendel, 28:23–31, 06 2022.
  * [19] M. Shoeybi, M. Patwary, R. Puri, P. LeGresley, J. Casper, and B. Catanzaro.  Megatron-lm: Training multi-billion parameter language models using model parallelism.  arXiv preprint arXiv:1909.08053, 2019.
  * [20] M. Tan and Q. Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In International conference on machine learning, pages 6105–6114. PMLR, 2019.
  * [21] M. Tan and Q. Le.  Efficientnetv2: Smaller models and faster training.  In International conference on machine learning, pages 10096–10106. PMLR, 2021.
  * [22] Z. Tu, H. Talebi, H. Zhang, F. Yang, P. Milanfar, A. Bovik, and Y. Li.  Maxvit: Multi-axis vision transformer.  In Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXIV, pages 459–479. Springer, 2022.
  * [23] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin.  Attention is all you need, 2017.
  * [24] Q. Wan, Z. Huang, J. Lu, G. Yu, and L. Zhang.  Seaformer: Squeeze-enhanced axial transformer for mobile semantic segmentation, 2023.
  * [25] C. Yu, C. Gao, J. Wang, G. Yu, C. Shen, and N. Sang.  Bisenet v2: Bilateral network with guided aggregation for real-time semantic segmentation.  International Journal of Computer Vision, 129:3051–3068, 2021.
  * [26] A. Zhang, Z. C. Lipton, M. Li, and A. J. Smola.  Dive into deep learning, 2023.
  * [27] Z. Zhang, A. Zhang, M. Li, H. Zhao, G. Karypis, and A. Smola.  Multimodal chain-of-thought reasoning in language models, 2023.
  * [28] S.-Y. Zhou and C.-Y. Su.  A novel lightweight convolutional neural network, exquisitenetv2, 2022\.