\section{Introduction}
Cardiovascular diseases (CDVs) are a group of disorders of the heart and blood vessels,  including coronary heart disease, cardiovascular disease, rheumatic heart disease and other conditions.
These diseases are the leading cause of global death. The
World Health Organisation (WHO) estimates that 32\% of all deaths worldwide are caused by CDVs. In addition to preventing CDVs by
maintaining a healthy lifestyle, early detection can curb mortality rates for heart disease patients and also for people with increased risk of this condition.

An electrocardiogram (ECG) is a safe and  non-invasive  procedure
  that  detects cardiac abnormalities by measuring the electrical activity generated by the heart as it contracts.
Sensors attached to the skin are used to detect the electrical signals produced by the heart each time it beats. These signals are
 analysed  by a  specialist to assess if they are unusual.
Computerised ECG analysis plays a critical
role in clinical ECG workflow. It is predicted  that in the
future full automated analysis of ECG may become more reliable
than human analysis.

Automation of ECG interpretation has been successfully
 addressed  by Artificial Intelligence (AI) methods
  which fall within the end-to-end deep learning
 framework \cite{AHO17,MR19, HRH19,RRP20}. However, due to the
over-parameterised black-box nature of this framework \cite{LBH15} it is difficult to understand how deep models make decisions.
Besides, interpretation
using Deep Neural Networks (DNN) has been shown to be
 susceptible to adversarial attacks \cite{HHL20,ML22}.
 Automatic classification of heartbeats
 has also been addressed by different techniques of
 machine learning which rely on the ability to
extract distinctive features of the beats in each class. Reviews and
summaries of these techniques can be found in
\cite{LSC16,AFA17,LMM18,MM22}.

In this work we present an alternative viewpoint for the
  specific task of
 morphological identification of heartbeats. We base the proposed
  method on the ability to model the morphology of
a heartbeat as superposition of elementary components. We work under the  assumption that there exits a set of elementary component which is
  better suited for representing heartbeats of certain morphology
 than others. Inspired by a previous  work on face recognition
 \cite{WYG08}, we consider that the sets of elementary components
are redundant. However, we differ from \cite{WYG08} in the methodology  and also
in that here the redundant sets are to be learned from  data.
 The common ground between the approaches is that both base
 the identification criterion on the concept of sparsity.

In the signal processing field, a signal approximation
is said to be sparse if it belongs to a
subspace of much smaller dimensionality than that of the space were the signal comes from. The search for the approximation subspace is carried out using  a redundant set, called a `dictionary'. Only a few
elements of the dictionary, called `atoms', are finally used for
constructing the signal approximation, called `atomic decomposition'.
 Sparse approximation of
heartbeats is applied in \cite{RR18} for classification using Gabor dictionaries.
The particular Gabor's atoms involved in the representation of
the beats, as well as the coefficients in the atomic decomposition,
are taken as features to feed machine learning classifiers.
Contrarily, our proposal relies on the possibility  of
learning dictionaries, which play themselves a central role in the decision
making process.  To the best of our knowledge this setup has not been considered before. We focus on the layout of the approach and
  test it for binary classification.

The paper is organised as follows. Sec.~\ref{SRUD} discuses the
greedy algorithms to be considered, as well as  the method
for learning the dictionaries. The sparsity metrics supporting the
discrimination criteria are also introduced in this section. The
 numerical tests illustrating the approach are
placed in Sec.~\ref{numex}.
The conclusions are drawn in
Sec.~\ref{conclu}.
\section{Signal representation using dictionaries}
\label{SRUD}
We start this section by introducing
some basic notation.
Throughout the paper
$\mathbb{R}$ represents the set of real numbers, and  $\mathbb{N}$ and
$\mathbb{Z}$ the sets of natural and integer numbers respectively.
Boldface lower case letters indicate Euclidean vectors and
boldface capital letters indicate matrices. The
corresponding components are indicated using standard mathematical fonts
e.g., $\mathbf{d} \in \mathbb{R}^N$ is a vector of components
$d(i)\in \mathbb{R},\, i=1,\ldots,N$ and $\mathbf{G} \in \mathbb{R}^{N\times M}$ is a matrix of components $G(i,j),\ i=1,\ldots, N, \, j=1,\ldots,M$.
The Euclidean inner product $\left \langle \cdot, \cdot \right \rangle$ between vectors
in $\mathbb{R}^N$ is defined as $$\left \langle \mathbf{d} ,\mathbf{g} \right \rangle= \sum_{i=1}^{N}
  d(i)g(i).$$
This definition induces the 2-norm
$\|\mathbf{g}\| = \sqrt{\left \langle \mathbf{g} ,\mathbf{g} \right \rangle}$. The 1-norm of a
vector $\mathbf{c} \in \mathbb{R}^{N}$, is indicated as $\|\mathbf{c}\|_1$ and
calculated as $\|\mathbf{c}\|_1= \sum_{n=1}^N |c(n)|$.

The inner product between matrices in $\mathbb{R}^{N \times M}$
is the Frobenius inner product, which is
defined as
$$\left \langle \mathbf{F} ,\mathbf{G} \right \rangle_F= \sum_{i=1}^{N} \sum_{j=1}^{M}
  F(i,j)G(i,j).$$
This definition induces the Frobenius norm
$\|\mathbf{G}\|_F = \sqrt{\left \langle \mathbf{G} ,\mathbf{G}\right \rangle_F}.$
The transpose of a matrix $\mathbf{G}$ is
indicated as $\mathbf{G}^\top$.

Given a signal $\mathbf{f}$ as a vector in $\mathbb{R}^N$,
the $K$-term
atomic decomposition for its
  approximation is
 of the form
\begin{equation}
\label{atoK}
\mathbf{f}^K=\sum_{n=1}^K c(n)\mathbf{d}_{\ell_n},
\end{equation}
where the elements $\mathbf{d}_{\ell_n}$, called atoms, are
chosen from a redundant set
$\mathcal{D}=\{\mathbf{d}_n, \in \mathbb{R}^N, \|\mathbf{d}_n\|=1\}_{n=1}^M$,
called dictionary. The set of all linear combination of
elements in $\mathcal{D}$ is denoted as ${\mbox{\rm{span}}}\{\mathcal{D}\}$.
The problem of how to select from $\mathcal{D}$ the
$K$ elements $\mathbf{d}_{\ell_n},\,n=1\ldots,K$ such that
$\|\mathbf{f}^K - \mathbf{f}\|$ is minimal is intractable (there
are $\frac{N!}{(N-K)!K!}$ possibilities).
In practical applications the problem is addressed by `tractable' methods. For the most part these methods are realised by
\begin{itemize}
\item [(a)] Expressing $\mathbf{f}^K=\sum_{n=1}^M c(n)\mathbf{d}_n$
 using only  $K$-nonzero coefficients minimising
 the 1-norm $\|\mathbf{c}\|_1=\sum_{n=1}^M  |c(n)|$ \cite{CDS98}.

\item [(b)] Using a greedy pursuit strategy for stepwise
selection of the $K$ elements $\mathbf{d}_{\ell_n},\,
n=1,\ldots,K$, for producing the approximation
\eqref{atoK}.
\end{itemize}

We restrict consideration to greedy pursuit algorithms,
because for the application we are considering these types of methods are effective and faster that those based on minimisation of the 1-norm.

\subsection{Pursuit Strategies}
In the context of signal processing the simplest pursuit strategy
 is known under the name of Matching Pursuit (MP)\cite{MZ92}.
Depending
on their
implementation and context of application
variations of
the  MP
approach can also be found under different names.
 We discuss here a refinement of MP known as
Orthogonal Matching Pursuit (OMP) \cite{PRK93} as well as
the stepped wise
optimised version termed
 Optimised Orthogonal Matching
 Pursuit (OOMP) \cite{RNL02}.

\subsubsection {From MP to OMP}
The MP algorithm evolves by successive approximations as follows:
 Setting $k=0$ and starting with an initial approximation
$\mathbf{f}^0=0$ and residual $\mathbf{r}^{0} = \mathbf{f}$,
the algorithm
  progresses  by sub-decomposing the $k$-th order residual
in the form
\begin{equation}
\mathbf{r}^{k} =
\left \langle \mathbf{d}_{\ell_{k+1}}, \mathbf{r}^{k} \right \rangle \mathbf{d}_{\ell_{k+1}}
+ \mathbf{r}^{k+1},
\label{tech:1}
\end{equation}
where $\mathbf{d}_{\ell_{k+1}}$ is the atom corresponding to the index selected as
\begin{equation}
\label{selMP}
\ell_{k+1}=  \operatorname*{arg\,max}_{\substack{n=1,\ldots,M}} |\left \langle \mathbf{d}_{n} , \mathbf{r}^{k} \right \rangle|.
\end{equation}
This atom is used to update the approximation $\mathbf{f}^k$ as
\begin{equation}
\label{upfk}
\mathbf{f}^{k+1} = \mathbf{f}^{k} + \left \langle \mathbf{d}_{\ell_{k+1}}, \mathbf{r}^{k} \right \rangle  \mathbf{d}_{\ell_{k+1}}.
\end{equation}
From \eqref{tech:1} it follows that
$\|\mathbf{r}^{k+1}\| \le \|\mathbf{r}^{k}\|$, since
\begin{equation}
\|\mathbf{r}^{k}\|^{2} = |\left \langle \mathbf{d}_{\ell_{k+1}}, \mathbf{r}^{k} \right \rangle|^{2} + \|\mathbf{r}^{k+1}\|^{2}.
\label{tech:2}
\end{equation}
It is easy to prove that in the limit $k \rightarrow \infty$, the sequence  $\mathbf{f}^{k}$
given in \eqref{upfk} converges to $\mathbf{f}$, if $\mathbf{f} \in V_{M}={\mbox{\rm{span}}}\{\mathcal{D}\}$, or to
$\hat{P}_{V_{M}}\mathbf{f}$ the orthogonal projection of $\mathbf{f}$ onto
$V_{M}$, if $\mathbf{f} \notin V_{M}$ (a pedagogical proof can be found on \cite{RNRS20}). However, the  method is not
 stepwise optimal because it does not
yield an orthogonal projection at each
step.  Accordingly, the algorithm may select linearly dependent atoms, which is the main drawback of MP when applied with highly coherent dictionaries.
A refinement to MP, which does yield an
orthogonal projection approximation at each step is called
OMP \cite{PRK93}. The method selects the atoms as in
\eqref{selMP} but at each iteration produces a decomposition of the signal as given by:
\begin{equation}
\mathbf{f}^k = \sum_{n=1}^{k} c^k(n) \mathbf{d}_{\ell_{n}} + \tilde{\mathbf{r}}^{k},
\label{tech:5}
\end{equation}
where the coefficients $c^k(n)$ are computed in such a way
that
it is true that
$$\sum_{n=1}^{k} c^k(n) \mathbf{d}_{\ell_{n}}= \hat{P}_{\mathtt{V}_{k}}\mathbf{f},\quad{\text{with}}
\quad V_{k}= {\mbox{\rm{span}}}\{\mathbf{d}_{\ell_{n}}\}_{n=1}^k.$$
The superscript of $c^k(n)$ in \eqref{tech:5} indicates
the dependence of these quantities on the iteration step $k$.
Thus,
in addition to selecting linearly independent atoms,
OMP yields the unique element $\mathbf{f}^k \in V_{k}$ minimising
$\|\mathbf{f}^k -\mathbf{f}\|$. As discussed  next,
a further refinement to OMP, called
OOMP \cite{RNL02}, selects also the atoms in
order to minimise in
a stepwise manner the norm of the residual error.
\subsubsection{The OOMP method}
The OOMP approach iterates as follows.
The algorithm is initialised by setting:
$\mathbf{r}^0=\mathbf{f}$, $ \mathbf{f}^0=0$, $\Gamma= \emptyset$
and $k=0$. The
 first atom
 is selected as the one corresponding to the index
$\ell_{1}$ such that
\begin{equation}
\ell_{1}=\operatorname*{\,max}_{n=1,\ldots,M}
 \left |\left \langle \mathbf{d}_n,\mathbf{r}^{0}\right \rangle \right|^2.
\end{equation}
This first atom is used to
assign $\mathbf{w}_1=\mathbf{d}_{\ell_{1}}=\boldsymbol{\beta}_1$,
calculate $\mathbf{r}^{1}= \mathbf{f} - \mathbf{d}_{\ell_{1}}
\left \langle \mathbf{d}_{\ell_{1}}, \mathbf{f}\right \rangle$, and  iterate
 as prescribed below.

\begin{itemize}
\item[1)] Upgrade the set $\Gamma \leftarrow  \Gamma \cup
\ell_{k+1}$, increase $k \leftarrow k +1$, and
 select the index of a new atom for the approximation as
\begin{equation}
\label{oomp}
\ell_{k+1}=\operatorname*{\,max}_{\substack{n=1,\ldots,M\\ n\notin \Gamma}}
 \frac{|\left \langle \mathbf{d}_n,\mathbf{r}^{k} \right \rangle|^2}{1 - \sum_{i=1}^{k}
|\left \langle \mathbf{d}_n ,\til{\mathbf{w}}_i\right \rangle|^2},
 \quad \text{with} \quad \til{\mathbf{w}}_i= \frac{\mathbf{w}_i}{\|\mathbf{w}_i\|_2}.
\end{equation}
\item[2)]
Compute the corresponding new vector $\mathbf{w}_{k+1}$ as
\begin{equation}
\begin{split}
\label{GS}
\mathbf{w}_{k+1}= \mathbf{d}_{\ell_{k+1}} - \sum_{n=1}^{k} \frac{\mathbf{w}_n}
{\|\mathbf{w}_n\|^2_2} \left \langle \mathbf{w}_n, \mathbf{d}_{\ell_{k+1}}\right \rangle.
\end{split}
\end{equation}
including, for numerical accuracy,  the
re-orthogonalisation step:
\begin{equation}
\label{RGS}
\mathbf{w}_{k+1} \leftarrow \mathbf{w}_{k+1}- \sum_{n=1}^{k} \frac{\mathbf{w}_{n}}{\|\mathbf{w}_n\|_2^2}
\left \langle \mathbf{w}_{n} , \mathbf{w}_{k+1}\right \rangle.
\end{equation}

\item[3)]
 Upgrade vectors $\boldsymbol{\beta}_n^{k}$ as
\begin{equation}
\boldsymbol{\beta}_{k+1}^{k+1}= \frac{\mathbf{w}_{k+1}}{\|\mathbf{w}_{k+1}\|^2}, \quad
\boldsymbol{\beta}_{n}^{k+1}= \boldsymbol{\beta}_{n}^{k} - \boldsymbol{\beta}_{k+1}^{k+1}\left \langle \mathbf{d}_{\ell_{k+1}}, \boldsymbol{\beta}_{n}^{k}\right \rangle,\quad n=1,\ldots,k.
\end{equation}
\item [4)] Upgrade vector $\mathbf{r}^{k}$ as
\begin{eqnarray}
\mathbf{r}^{k+1} &=& \mathbf{r}^{k} - \left \langle \mathbf{w}_{k+1}, \mathbf{f} \right \rangle  \frac{\mathbf{w}_{k+1}}{\|\mathbf{w}_{k+1}\|^2}.
\end{eqnarray}
\item[5)] If for a given $\rho$ value
  the condition
 $\|\mathbf{r}^{k+1}\| < \rho$ has been met
 stop the selection process. Otherwise repeat steps 1) - 5).
\end{itemize}
Calculate the coefficients  in \eqref{atoK} as
$$\mathbf{c}=\left \langle \boldsymbol{\beta},\mathbf{f}\right \rangle.$$
Calculate the final approximation $\mathbf{f}^k$ as
$$\mathbf{f}^k = \mathbf{f} -\mathbf{r}^k.$$
\begin{remark}
The difference between  our implementations of OMP and
OOMP is only that for OMP the denominator in
\eqref{oomp} is eliminated.  The
selection is realised as in MP (c.f.  \eqref{selMP})
which is not a stepwise optimal selection. Instead,
	\eqref{oomp}  selects the atom $\mathbf{d}_{\ell_{k+1}}$ which, fixing the previously selected atoms $\mathbf{d}_{\ell_i},\,i=1,\ldots,k$,
	minimises the distance $\|\mathbf{f}^{k+1} -\mathbf{f}\|$. The proof can be found in \cite{RNL02}.
\end{remark}
\subsection{Sparsity criterion for morphological feature
extraction of heartbeats}
\label{SM}
A digital ECG signal is a
 sequence of
heartbeats, each of which is characterised by a
combination of three graphical deflections, known
as QRS complex, and the so called P and T waves.
An idealised illustration of these deflections are
given in Fig.~\ref{qrs}.

\begin{figure}[h]
\centering
        \includegraphics[height=.4\textheight]{SinusRhythmLabels.eps}
\caption{Public domain graph of the QRS complex and P and T waves representing a normal heartbeat.}
\label{qrs}
\end{figure}
However, as illustrated in Fig.~\ref{beats}, in real ECG
records the shape of the beats in the same class may vary.
\begin{figure} [!h]
\centering
        \includegraphics[height=.2\textheight]{Beat_N1.eps}
        \includegraphics[height=.2\textheight]{Beat_N2.eps}
\caption{Two normal heartbeats from different
ECG records}
\label{beats}
\end{figure}

Approximation for morphological feature extraction of a heartbeat
requires the segmentation of the QRS complex.
In this work the segmentation is realised
by taking  a fixed number of samples to left and right of
the location of the  R-peak, to have beats
of equal support, say $N_q$.

Assuming that one wishes to morphologically differentiate heartbeats of, say class $\mathtt{A}$, from heartbeats of class $\mathtt{B}$, and that
a dictionary $\mathcal{D}_{\mathrm{a}}$, specially
constructed to approximate a heartbeat of
class $\mathtt{A}$, and a dictionary $\mathcal{D}_{\mathrm{b}}$,  specially
constructed to approximate a heartbeat
of class  $\mathtt{B}$, are given, we discuss next
several possibilities of using sparsity as a distinguishability
criterion.

In order to decide whether a heartbeat
 $\mathbf{f} \in \mathbb{R}^{N_q}$ belongs to class $\mathtt{A}$ or class $\mathtt{B}$ this
beast is approximated, up to the same precision, using
both dictionaries, i.e.
\begin{equation}
\label{K1}
\mathbf{f}^{a}=\sum_{n=1}^{K_a} c^a(n)\mathbf{d}_{\ell_n^a},
\end{equation}
where the atoms
 $\mathbf{d}_{\ell_n^a}$ are
chosen from $\mathcal{D}_{\mathrm{a}}$ and
\begin{equation}
\label{K2}
\mathbf{f}^{b}=\sum_{n=1}^{K_b} c^b(n)\mathbf{d}_{\ell_n^b},
\end{equation}
where the atoms $\mathbf{d}_{\ell_n^b}$ are
chosen from $\mathcal{D}_{\mathrm{b}}$. \\
\\
\subsubsection*{Criterion I (a) and (b) }
\label{crit}
If $K_a < K_b$ the beat $\mathbf{f}_q$ is assigned to class $\mathtt{A}$.\\
If $K_b < K_a$ the beat $\mathbf{f}_q$ is assigned to class $\mathtt{B}$.\\
If $K_a=K_b$ the criterion does not make a decision.\\

In the event that $K_a=K_b$ a decision could still be made
by recourse to a different metric.\\

\noindent
(a) {\em{Smaller entropy criterion:}}
Assigning $p^a(n)=\frac{|c^a(n)|}{\|\mathbf{c}^a\|_1}$ and $p^b(n)=\frac{|c^b(n)|}{\|\mathbf{c}^b\|_1}$ calculate the corresponding Shannon's entropies
$$S^a= -\sum_{n=1}^{K_a} p^a(n) \ln(p^a(n)),$$
$$S^b= -\sum_{n=1}^{K_b} p^b(n) \ln(p^b(n)).$$
The entropy of, say $p^a(n),\,n=1,\ldots,K_a$, would be smaller
if the components are fewer and the magnitude of
some  $p^a(n)$ much larger than others.
In particular the minimum entropy value occurs
if $p^a(n)=\frac{|c^a(n)|}{\|\mathbf{c}^a\|_1}=1$ for some
$n$ and zero otherwise, in which case $S^a=0$. Accordingly, if
$S_a < S_b$ the beat $\mathbf{f}$ might be assigned to class $\mathtt{A}$
and if $S_b < S_a$ the beat $\mathbf{f}$ might be assigned to class $\mathtt{B}$.
Let us recall that the criterion of smaller entropy was
 introduced for basis selection in
the context of wavelet packets \cite{CW92}.\\

\noindent
(b) {\em{Smaller norm-1 criterion.}}\\
This criterion is in line with the
{\em{basis pursuit approach}} \cite{CDS98}
which adopts the minimisation of the norm-1 as a way of
producing a tractable sparse solution from a given
dictionary.  Consequently, if
 $\|\mathbf{c}^a\|_1 < \|\mathbf{c}^b\|_1$ the beat $\mathbf{f}$ is assigned to
class $\mathtt{A}$
and if $\|\mathbf{c}^b\|_1 < \|\mathbf{c}^a\|_1$ the beat $\mathbf{f}$ if be assigned to class $\mathtt{B}$.
\subsubsection*{Criterion II}
Use always the smaller entropy criterion.
\subsubsection*{Criterion III}
Use always the smaller norm-1 criterion.
\subsection{Dictionary learning}
So far we have assumed to know the dictionaries
$\mathcal{D}_{\mathrm{a}}$ and $\mathcal{D}_{\mathrm{b}}$, best suited for each class of heartbeat.
In this section we discuss how these dictionaries can be learnt from a data set of annotated heartbeats.

The adopted  strategy to learn a dictionary, as a matrix $\mathbf{D}$,
using the heartbeats placed in an array $\mathbf{F} \in \mathbb{R}^{N_q \times Q}$
 proceed
through a 2 step process.\\

\noindent
{\em{Step 1}}
\begin{itemize}

\item Set $k=0$. Given an initial dictionary,  as a
matrix $\mathbf{D}^k \in \mathbb{R}^{N_q \times M}$, use a greedy algorithm to approximate
 each beat  $\mathbf{f}_q \in \mathbb{R}^{N_q}$, for  $q=1,\ldots,Q$,
as in \eqref{K1} or \eqref{K2}  and collect the  vectors of
 coefficients $\mathbf{c}_q^k \in \mathbb{R}^{K_q}$.

\item Place each
	vector $\mathbf{c}_q^k$ as a column of a matrix $\mathbf{C}^k \in \mathbb{R}^{M\times Q}$ having nonzero elements $C{({\ell_n^q},q)}^k={c_q(n)}^k,\,n=1,\ldots,{K_q}^k,\,
q=1,\ldots,Q$ and  set
 all the other elements equal to zero.

\item Using matrix $\mathbf{C}^k$ and dictionary $\mathbf{D}^k$ calculate
the approximation $\til{\mathbf{F}}$ of $\mathbf{F}$ as
$$\til{\mathbf{F}}^k =\mathbf{D}^k  \mathbf{C}^k$$

\end{itemize}

\noindent
{\em{Step 2}}
\begin{itemize}

\item Using the approximation $\til{\mathbf{F}}^k$ find the updated dictionary $\mathbf{D}^{k+1}$ minimising
	$\|\til{\mathbf{F}} - \til{\mathbf{F}}^k\|_F^2$.
Thus,
\begin{equation}
\label{dic}
\mathbf{D}^{k+1}=
		\argmin_{\substack{\mathbf{D}}}\|{\til{\mathbf{F}}} -  {\mathbf{D}} \mathbf{C}^k \|_F^2 = {\til{\mathbf{F}}} {\mathbf{C}^k}^\top (\mathbf{C}^k{\mathbf{C}^k}^\top)^{-1}.
\end{equation}

\item Given a maximum  number of iterations,  $\Max$ say,
and a tolerance $\mathrm{tol}$  for the error norm, if  $\|\mathbf{D}^{k+1} - \mathbf{D}^{k}\|_F < \mathrm{tol}$  or  $k+1 >\Max$  stop. Otherwise set $k \to k+1$ and repeat {\em{Steps 1}} and {\em{2}}.
\end{itemize}

\begin{remark}
The dictionary updating equation \eqref{dic}  requires
	that matrix $\mathbf{C}^k{\mathbf{C}^k}^\top$ should have an inverse. This would not be true if some elements in the dictionary were not chosen in the previous step. In that case, the unselected atoms
should be
removed from the dictionary before implementing equation
\eqref{dic}. In practice, as long as
$Q$ is significantly larger than $M$ and the examples are
independent,
	matrix $\mathbf{C}^k{\mathbf{C}^k}^\top$ has an inverse.
\end{remark}

\begin{remark}
 The problems of determining matrix $\mathbf{C}^k$ at Step 1
and  matrix $\mathbf{D}$ at Step 2 are
convex problems with unique solution. However, the
combined problem of determining matrix $\mathbf{C}^k$ and $\mathbf{D}$
 is not jointly convex. Thus, the solution depends on
the initial dictionary. As will be demonstrated by the
 simulations, this in not crucial within the context of
 the proposed approach.
\end{remark}

\section{Binary morphological differentiation of
heartbeats}
\label{numex}
In this section we test the proposal by differentiating the
classes ${\tt{N}}$ (Normal) and ${\tt{V}}$ (Ventricular Ectopic)
in the MIT-BIH Arrhythmia data set \cite{MITDB, MM01}.
As per the recommendations given by the Advancement of Medical Instrumentation (AAMI) the MIT-BIH Arrhythmia database is projected into five classes. Within these classes the ${\tt{N}}$ beats considered here
include: The normal beats N, the
Left Bundle Branch Block Beats (LBB)   and
the Right  Bundle Branch Block  Beats (RBB).
The ${\tt{V}}$ class comprises: the Premature Ventricular Contraction (PVC)  and Ventricular Escape Beats (VEB).
While the number of ${\tt{V}}$ beats is much less
than the number of ${\tt{N}}$ beats, the former are
still enough to learn the corresponding dictionary.
The locations of the R-peaks are retrieved from the
annotations provided with the dataset, using
the Matlab software \cite{GAG00} available on \cite{MITDB}. The peaks are then segmented by taking $145$ samples to the right and
110 samples to  the left of the
 R-peak location.
 Before learning the dictionary a quick check of the
 training set is realised,  in order to find out
 if there are some peaks
 that would not qualify as
 `good' examples. For this end we proceed as explained below.
\subsubsection*{Screening the data training sets}
For automatic screening of a set
the segmented beats of the same class
are approximated, up to the same quality,
using a greedy algorithm and a general wavelet dictionary.
We illustrate the process by giving the details for screening
the training set for the numerical Test I.

Wavelet dictionaries arise from translations,  by
a parameter $b =2^{-l},\,l \in \mathbb{N}$  of
scaling prototypes  \cite{ARN08,RNC19}
\begin{equation} \label{scaling_functions_dictionary}
\phi_{j_0,k,b} \left( x \right)  =  \phi  \left( 2^{j_0} x - b k \right), \quad k \in \mathbb{Z},
\end{equation}
and wavelet prototypes at different scales
\begin{equation} \label{wavelets_dictionary}
 \psi_{j,k,b} \left( x \right) =  \psi \left( 2^j x - b k \right), \quad
k \in \mathbb{Z}, \quad j \geq j_0.
\end{equation}
For screening of the training set we have used
the 9$/$7 Cohen-Daubechies-Feauveau ({\bf{cbf97}}) wavelet
family with $b=0.25$, $j_0=2$ and $j=3,4,5,6$,
which introduces a redundancy factor  of 2.67.
 This dictionary was
generated with the software described in \cite{CRN21} available on
\cite{webpage1}. The scaling and wavelet prototypes for
 the {\bf{cbf97}} wavelet family
 are shown in Fig.~\ref{cdf}.

\begin{figure} [!ht]
\centering
        \includegraphics[height=.2\textheight]{phi.eps}
        \includegraphics[height=.2\textheight]{psi.eps}
\caption{Graph of the scaling prototype $\phi(x)$ (left)
and wavelet prototype $\psi(x)$ (right)
	for the {\bf{cbf97}} family.}
  \label{cdf}
\end{figure}

The quality in the approximation
of each beat is fixed
using the ${\mathrm{prdn}}$ metric as defined  by
$${\mathrm{prdn}}(q)= \frac{\|\mathbf{f}_q - \mathbf{f}_q^k\|}{\|\mathbf{f}_q - \ov{\mathbf{f}}_q \|} \times 100\%, \quad q=1,\ldots,Q,$$
where each $\mathbf{f}_q$ is a heartbeat, $\mathbf{f}_q^k$ its corresponding approximation by $k$ atoms, and $\ov{\mathbf{f}}_q$ the mean value of
$\mathbf{f}_q$.
The approximation of different beasts, up to the same ${\mathrm{prdn}}$, is achieved for different values of $k$.
The left graph of Fig.\ref{his1} shows the histograms of
$k$ values obtained when approximating, up to
${\mathrm{prdn}}=9$, the ${\tt{N}}$ beats in the training set of the numerical
Test I (c.f. Sec.~\ref{TestI}).
\begin{figure} [!ht]
\centering
        \includegraphics[width=0.5\textwidth,height=.25\textheight]{HistN2.eps}
        \includegraphics[height=.25\textheight]{outlyingN2.eps}
  \caption{Histogram of $k$ values for approximating
	(up to ${\mathrm{prdn}}=9$)
 the ${\tt{N}}$-signals
in the training set (left graph).
Signal corresponding to the outlying value $k=63$ (right graph).}
  \label{his1}
\end{figure}

Even if not very noticeable in the histogram there are
a few values of $k$ very far from the main support.
These values are from `rare' signals that should be
 investigated. In this set they are just signals looking as pure
 noise, e.g. the signal in the right graph of Fig \ref{his1}
corresponds to
$k=61$, very far from the mean  value, $\overline{k}=14.31$
with standard deviation ${\mathrm{std}}=5.90$.
When learning the dictionary for
classification of ${\tt{N}}$ beats we disregard beats in the
training set producing values of $k$ outside the range
 $(\overline{k} - 3\, {\mathrm{std}}, \overline{k} + 3\, {\mathrm{std}})$.
This amounts to disregarding $2.58\%$
of the total ${\tt{N}}$ beats in the training set set.
\begin{figure} [h!]
\centering
        \includegraphics[width=0.47\textwidth,height=.25\textheight]{HistV.eps}
        \includegraphics[height=.25\textheight]{outlyingV2.eps}
  \caption{Histogram of $k$ values for  approximating
	(up to ${\mathrm{prdn}}=9$) the ${\tt{V}}$-signals
in the training set (left graph).
Signal corresponding the outlying value $k=41$ (right graph).
}
  \label{his2}
\end{figure}

The left graph of Fig.\ref{his1} shows the histograms of
$k$ values obtained when approximating, up to
${\mathrm{prdn}}=9$, the ${\tt{V}}$ beats in the training set of the numerical Test I.
For this class  $\overline{k}=10.61$ with ${\mathrm{std}}= 3.28$.
 Because the  data set of ${\tt{V}}$ beats is
 smaller than the previous one,
when learning the dictionary for
classification of ${\tt{V}}$ beats we disregard beats in the
training set outside the range
 $(\overline{k} - 2\,{\mathrm{std}}, \overline{k} + 2\,{\mathrm{std}})$.
This amounts to disregarding $2.65\%$
of the ${\tt{V}}$ beats in the training set.

\subsection{Numerical Test I}
\label{TestI}
The purpose of this test is to assess the suitability  of
binary morphological differentiation of
heartbeats using dictionaries learned from
examples of ${\tt{N}}$ and ${\tt{V}}$ shapes
in the dataset. For this end we randomly split the
 $Q=89840$  ${\tt{N}}$ beats into two groups:
 35\% of the ${\tt{N}}$ beats,
are used for training the dictionary, $\mathbf{D}_n$ say.
The remaining 65\% of the ${\tt{N}}$ beats are reserved
for testing.
Since the ${\tt{V}}$ beats are much less than the
${\tt{N}}$ ones, $50\%$ of the ${\tt{V}}$ beats are
used for training the dictionary $\mathbf{D}_v$, and the remaining
 $50\%$ for testing (c.f. Table~\ref{Table1})

\begin{table}
\begin{center}
\begin{tabular}{|l|c |c || c| c||}
\hline
Sets & \multicolumn{2}{|c||}{\rm{Training}} &
\multicolumn{2}{|c||}{\rm{Testing}}\\ \hline \hline
Class  & ${\tt{N}}$ & ${\tt{V}}$ & ${\tt{N}}$ & ${\tt{V}}$ \\ \hline \hline
        Number of beats & 30000 & 3359 & 57277&  3359\\ \hline \hline
\end{tabular}
\caption{Number of heartbeats of classes ${\tt{N}}$ and ${\tt{V}}$
in the training and testing sets for Test I.}
\end{center}
\label{Table1}
\end{table}

The dictionaries are learned
to have redundancy two,
i.e. each dictionary is a matrix of real numbers of
 size $256 \times 512$. We test the method against:\\
(i) The initial dictionary.\\
(ii) The different greedy algorithms considered in this
work: MP, OMP, OOMP.\\

For this test
we randomly take 512  beats from the $30000$ ${\tt{N}}$
 heartbeats in the training set to construct a matrix $\mathbf{D}_n$
of size
 $256 \times 512$, which is used as initial dictionary.
The learning curves for the dictionary $\mathbf{D}_n$
with the 3 greedy algorithms are shown in left graph of Fig.\ref{lcn}.
We repeat the process but taking the 512
dictionary atoms randomly from the
$3358$ ${\tt{V}}$ beats in the training set. The
learning curves for this dictionary are
  shown in the right graph of Fig.\ref{lcn}.

 \begin{figure} [!ht]
\centering
        \includegraphics[height=.25\textheight]{L_curveNT.eps}
        \includegraphics[height=.25\textheight]{L_curveVT.eps}
  \caption{Dictionary learning curves with the tree greedy algorithms.
	  The left graph corresponds to dictionary $\mathbf{D}_n$ and
	  the right graph to dictionary $\mathbf{D}_v$.}
  \label{lcn}
\end{figure}

The classification performance is assessed by  means of the
 true positive (TP), false positive (FP), and false negative (FN)
 outcomes. These values
are used to calculate the following statistics metrics
for each class.\\

Sensitivity ($\mathrm{SE}$): Number of correctly classified heartbeats
among the total number of beats in the set, i.e.

$$\mathrm{SE}= \frac{\mathrm{TP}}{\mathrm{TP} + \mathrm{FN}}.$$

Positive predictivity (PP): Ratio of correctly classified heartbeats   to all the beats classified in that class.

$$\mathrm{PP}= \frac{\mathrm{TP}}{\mathrm{TP} + \mathrm{FP}}.$$

Additionally, the total accuracy ($\mathrm{AC}$) of the classification is
calculated as the fraction of correctly classified heartbeats in
both classes

$$\mathrm{AC}=\frac{\text{Correctly classified beats in both classes}}{\text{Total number of beats}}.$$

In Table \ref{StT} these scores are given
 as the mean value
of 5 realisations corresponding to 5 random initialisations
in the dictionary learning  process. The
standard deviations (std) are shown in rows 4,6,8,10, and 12
of Table \ref{StT}. The Matlab codes for implementing the approach
with the all the tree greedy algorithms have been made available
on \cite{webpage2}.\\

\begin{table}[!ht]
\begin{center}
\begin{tabular}{|l||c|c|c|c||c|c|c|c||c|c|c|c||}
\hline

	& \multicolumn{4} {|c||}{MP}& \multicolumn{4} {c||}{OMP}&
	\multicolumn{4} {|c||}{OOMP} \\\hline \hline
	Crit.	& I (a)& I (b) & II & III& I (a) & I (b) & II & III&
	 I (a)& I (b) & II & III\\ \hline \hline
	$\mathrm{SE}_{{\tt{N}}}$ (\%)&{\bf{99.7}}&{\bf{99.7}}&{\bf{99.7}}& {\bf{99.7}}&99.1&99.1 &99.3&{\bf{99.4}}&{\bf{98.8}}&{\bf{98.8}}&98.5& 97.4\\\hline
std &0.05&0.05&0.04&0.03&0.09&0.09&0.01&0.04&0.14&0.13&0.16&0.45\\ \hline \hline
	$\mathrm{SE}_{{\tt{V}}}$ (\%)&95.7&95.7&{\bf{96.9}} &{\bf{97.6}}&
	96.5&96.5&96.7&{\bf{97.0}}&
	{\bf{97.6}}&{\bf{97.6}}&97.3&95.3\\\hline
        std &0.65 &0.63&0.45 &0.34&
	0.24&0.18&0.27&0.25&
	0.35&0.21&0.47&1.17\\ \hline \hline
	$\mathrm{PP}_{{\tt{N}}}$ (\%)&99.7&99.7&99.8&{\bf{99.9}}& {\bf{99.8}} &{\bf{99.8}}&{\bf{99.8}}&{\bf{99.8}}&{\bf{99.9}}&{\bf{99.9}}&99.8&99.7\\ \hline
        std &0.04&0.04&0.03&0.02&0.01&0.01 &0.02&0.01&0.02&0.01&0.03&0.07\\ \hline \hline
	$\mathrm{PP}_{{\tt{V}}}$ (\%)&{\bf{94.9}}&{\bf{94.9}}&94.6&94.1&86.4&86.6&{{88.6}}&{\bf{90.4}}&{\bf{82.5}}&{\bf{82.5}}&79.5&68.4\\\hline
        std &0.71&0.71&0.55&0.48&1.27&1.27&0.17&0.55&1.56&1.72&1.72&3.48\\ \hline \hline
	$\mathrm{AC}$ (\%)&{\bf{99.5}}&{\bf{99.5}}&{\bf{99.5}}&{\bf{99.5}}&99.0&{99.0}&99.1&{\bf{99.3}}&{\bf{98.7}}&{\bf{98.7}}&98.5&97.3\\\hline
        std &0.03&0.03&0.02&0.03&0.01&0.01&0.02&0.05&0.12&0.13&0.12&0.38\\ \hline \hline
\end{tabular}
	\caption{Average statistics scores (of 5 different initialisations in
	the dictionary learning procedure)
	 for classes ${\tt{N}}$ and ${\tt{V}}$, obtained with MP, OMP, and OOMP,  and the  sparsity criteria discussed in Sec.~\ref{SM}. std indicates the corresponding standard deviations and  $\mathrm{AC}$ the average accuracy.}
	 \label{StT}
	\end{center}
\end{table}

\newpage
{\bf{Discussion of results}}\\
\noindent

\noindent
i) The high accuracy and $\mathrm{SE}$ values for both classes,
indicates that:
35\% of the whole ${\tt{N}}$  beats in the data set and
50\% of the whole ${\tt{V}}$ beats in the data set provide enough
examples in the training set to learn dedicated dictionaries for each class. \\

\noindent
ii) From Table I we can assert that each of the greedy algorithms
performs better with a particular decision criterion.\\

\noindent
iii) MP, combined with the decision Criterion II (norm-1), produces
the highest statistics scores. \\

\noindent
iv) For most of the scores, the random initialisation of the learning process does not change the results in any significant manner (low std values). \\

\noindent
v) The comparatively much lower values  of $\mathrm{PP}_{{\tt{V}}}$
 reflex the enormous difference of the
number of ${\tt{N}}$  and ${\tt{V}}$ beats in the testing sets
(57277 vs  3359  beats). Hence, even if the percentage of incorrectly classified ${\tt{N}}$ beats is small, the classification being
binary implies that the number of incorrectly classified ${\tt{N}}$ beats count as false positive ${\tt{V}}$. This causes the low $\mathrm{PP}_{{\tt{V}}}$ values in comparison to the other  statistics metrics.
\\

\noindent
vi) Even if the OOMP approach secures the least number of atoms in th learning process (c.f. Fig.\ref{lcn}) it is not the approach rendering the highest classification scores. This is because  what matters for classification is the relative sparsity with respect to the
2 dictionaries. It is not surprising, though, that OOMP is the only approach that works better with Criterion I.

\begin{remark}
Accuracy of 99\% for differentiating classes in the whole
	MIT-BIH Arrhythmia  data set is state of the art result obtainable from
	other features and  other
	 machine learning classifiers. See for instance\\
	 \begin{itemize}
	 \item[\cite{MAP13}] $\mathrm{AC}\,\, 99\%$for classification
			 of 2 classes,
	 \item[\cite{AS15}] $\mathrm{AC} \,\,99\%$  for classification  of 5 classes,
	 \item[\cite{RRS16}] $\mathrm{AC}\,\, 99.2\%$ for classification  of 16 classes,
	 \item[\cite{RR18}] $\mathrm{AC}\,\, 99.3\%$ for classification  of 16 classes.
	 \end{itemize}
Although results are not strictly comparable, because the number of
classes  are different, we believe it appropriate to highlight
	that the
99.5\% accuracy attained in this work achieves equivalent
	values as other approaches do for binary and
	multi classification.
\end{remark}
\noindent
\subsection{Test II}
We test now the proposal in a more challenging situation.
Because the morphology of ${\tt{N}}$ and ${\tt{V}}$ depends on the
particular ECG signal, a realistic test is carried out
 by taking the training set and testing set from different records (corresponding to different patients).

Following previous publications \cite{COR04, RR18}
the data are divided for training and testing as below.

The training set is taken only from the 21 records below
in the the MIT-BIH Arrhythmia  data set.

\noindent
101, 106,  108, 112, 114, 115, 118, 119, 122, 124, 201, 203, 205, 207, 208, 209, 215, 220, 223, 230.\\

The remaining 21 records provide the testing set. These are
the records\\

\noindent
100, 103, 105, 111, 113, 117, 121, 123, 200, 202, 210, 212, 213, 214, 219, 221, 222, 231, 232, 233, 234
in the  MIT-BIH Arrhythmia  data set.\\

\begin{table}[!ht]
\begin{center}
\begin{tabular}{|l||c|c|c|c||c|c|c|c||c|c|c|c||}
\hline
        & \multicolumn{4} {|c||}{MP}& \multicolumn{4} {c||}{OMP}&
        \multicolumn{4} {|c||}{OOMP} \\\hline \hline
        Crit.   & I (a)& I (b) & II & III& I (a) & I (b) & II & III&
         I (a)& I (b) & II & III\\ \hline \hline
	$\mathrm{SE}_{{\tt{N}}}$ (\%)&89.8 & 89.9  & 90.5 & {\bf{91.8}}
	&88.8 & 89.5  & 88.8 & {\bf{92.5}}
	&88.7 & {\bf{89.0}}  & 85.7 & 81.5\\ \hline
std	&1.43 &  1.39 & 1.38 &1.26
	&1.39 &  1.28 & 1.48 &1.19
	&0.92 &  0.87 & 0.98 &2.10 \\ \hline \hline
$\mathrm{SE}_{{\tt{V}}}$ (\%) &87.6 & 87.6 & 90.2 & {\bf{91.0}}
	&89.8  &89.5 & {\bf{90.5}} & 88.9
	&92.3 & {\bf{92.6}} & 89.6 & 84.7 \\ \hline
std	&1.57& 1.64& 0.75& 1.14
	&1.40& 1.24& 1.46& 1.14
	&0.76& 0.94& 1.15& 2.83 \\ \hline \hline
	$\mathrm{PP}_{{\tt{N}}}$ (\%)&99.0& 99.0& 99.2& {\bf{99.3}}
	&{\bf{99.2}}&{\bf{99.2}}& {\bf{99.2}}& 99.1
	&{\bf{99.4}}& {\bf{99.4}}& 99.1& 98.7 \\ \hline
std	&0.11& 0.12& 0.05& 0.08
	&0.11& 0.09& 0.11& 0.08
	&0.06& 0.07& 0.10& 0.22 \\ \hline \hline
$\mathrm{PP}_{{\tt{V}}}$	&38.8 & 39.1& 41.1 &{\bf{45.1}}
	&37.1 & 38.5& 37.3 &{\bf{46.6}}
	&37.5 & {\bf{38.2}}& 31.5 &25.2 \\ \hline
std	&3.01& 2.96 &3.30 &3.56
	&2.70& 2.70 &2.73 &4.10
	&1.96& 1.85 &1.47 &1.69 \\ \hline \hline
$\mathrm{AC}$	&89.7& 89.8 &90.4 &{\bf{91.8}}
	&88.9& 89.5 &88.9 &{\bf{92.2}}
	&89.0& {\bf{89.3}} &86.0 &81.7 \\ \hline
std	&1.25& 1.21 &1.24 &1.10
	&1.24& 1.14 &1.33 &1.08
	&0.84& 0.79 &0.91 &1.81 \\ \hline \hline
\end{tabular}
\caption{Same description as in Table~\ref{StT} but for the
	numerical Test II, in which
training and testing sets correspond to different
	records.}
\label{StT2}
\end{center}
\end{table}

\noindent
{\bf{Discussion of results}}\\

\noindent
i) As expected, the classifications scores are lower than in
Test I. This is because not all the shapes of the beats in the
testing set
bears similarity with shapes present in the training set. \\

\noindent
ii) Also in this test each of the greedy algorithms
performs better with a particular decision criterion.\\

\noindent
iii) The greedy algorithm OOMP is the only one that performs
better with Criterion I.\\

\noindent
iv) Clearly MP works best with Criterion III  and for  most scores
OMP as well. \\

\noindent
v) The low values of $\mathrm{PP}_{{\tt{V}}}$ have the same cause as in Test I.
Here all values are even lower because the incorrectly classified
${\tt{N}}$ beats are more than in Test I.\\

\begin{remark}
In order to put the results into context we give here some
	indication of
	other technique performance. Nonetheless, it should be stressed once again that scores are not strictly comparable due to
	difference in the number of classes and, for the binary case,
	the classes and records for learning and testing.
\begin{itemize}
	\item[\cite{AY19}] $\mathrm{AC} \,\, 97\%$ for binary classification with 37 records for
		training and 11 records for testing, two of which are from a different data set.
	\item[\cite{ZLC14}] $\mathrm{AC}  \,\,86.7\%$ for classification of 4 classes.
	\item[\cite{RRS16}] $\mathrm{AC}  \,\,89.18 \%$ for classification of 5 classes.
	\item[\cite{RR18}]  $\mathrm{AC}  \,\,89.93 \%$ classification in of 5 classes.
\end{itemize}
\end{remark}
\section{Conclusions}
\label{conclu}
A set up for binary morphological identification of
heartbeats on the basis of sparsity metrics has been laid out.
The proposal was tested for identification of ${\tt{N}}$ and ${\tt{V}}$
beats in the MIT-BIH Arrhythmia data set, achieving
state of the art scores.

Because the number of ${\tt{N}}$ and ${\tt{V}}$ beats in the
 MIT-BIH Arrhythmia are unbalanced,
 the $\mathrm{AC}$ score is not reliable for assessment. Nevertheless,
the percentage of correctly classified ${\tt{N}}$ and ${\tt{V}}$ are
similar to the $\mathrm{AC}$ score
 in both numerical tests.

The results are encouraging because the proposed binary
identification is realised outside the usual machine leaning
framework, using sparsity as a
{\em{single}} parameter for making a decision. Thus, extensions of
the approach to allow for combination with other features and
other  machine leaning techniques are readily foreseen.

 The possibility of implementing the technique depends on the
 availability of `enough' examples  to learn the
 corresponding dictionaries. In the numerical tests
 realised here, for instance, 3359 examples for beats ${\tt{V}}$ were
 enough to learn the dictionary $\mathbf{D}_v \in \mathbb{R}^{256 \times 512}$.


\section{Introduction}\label{Introduction}

The purpose of this paper is to examine the meso-scale limit of a discrete micro-level mathematical model constructed to represent simple heat transfer.  The creation of the microscopic model is based on the application of the principle of dynamic consistency (\cite{13,  14,  15}).  Under the appropriate mathematical assumptions, we are able to obtain the previous results of  \cite{11}, \cite{4}, and \cite{16} regarding the replacement of the standard heat equation
\begin{equation}
u_t=Du_{xx},\label{eq1}
\end{equation} by the generalization
\begin{equation}
\tau u_{tt}+u_t=Du_{xx},  \label{eq2}
\end{equation}
where $D$ is the temperature diffusion constant,  $\tau$ is a time-lag parameter,  and
\begin{equation*}
u_t = \frac{\partial u}{\partial t},  \quad u_{tt} = \frac{\partial^2u}{\partial t^2},  \quad u_x=\frac{\partial u}{\partial x},  \quad u_{xx}=\frac{\partial^2 u}{\partial x^2}.
\end{equation*}
Note that $(x,t)$ are the one-dimensional space and time independent variables, and $u=u(x,t)$ is defined over appropriate intervals of $(x,t)$ with suitable boundary conditions and initial values. The need for a generalization of Eq. ~(\ref{eq1}) comes from the fact that the solutions of Eq.~(\ref{eq1}) transmit information at an infinite speed (\cite{2,5,6,7,8}),  a condition which violates the principle of causality (\cite{2,5,6,7,8}).

For convenience,  we now provide a resume of our major results:
\begin{itemize}
\item[(i)] Based on our formulation of the discrete mirco-level model, the drived first modified form of the standard heat PDE has the mathematical structure
\begin{equation} \label{eq3}
\epsilon_1 u_{tt} +u_t= Du_{xx}+(\epsilon_2D)u_{xxxx}.
\end{equation}
where $(\epsilon_1, \epsilon_2)$ are non-negative parameters.
\item[(ii)] Finite-order approximate models of discrete micro-level models can never be ``exact." By this, we mean that the meso-level equations will always require knowledge of certain properties of the solutions which in practice and/or theoretically can not be experimentally measured.
\end{itemize}

This paper is organized as follows: Section 2 provides a brief historical summary of work on the standard heat PDE,  its major defect,  and previous attempts to resolve the associated issue. It also includes several of the mathematical models put forth to solve the causality problem and discusses why they all have failed in this quest. Section 3 gives several general comments on the modeling process and introduces the concept of dynamic consistency.  In Section 4, we use dynamic consistency to construct perhaps the simplest discrete mathematical model for heat conduction; it corresponds to a generalized random walk model.  We take this to be our micro-level model. In Section 5, we show that by using a Taylor series expansion, a meso-level hybrid difference-differential equation is obtained.  Special cases of this equation includes all previous cases of the modified standard heat equation.  Finally,  in Section 6,  we summarize our major conclusions and provide reasons why the usual methods for generalizing the heat PDE may never produce valid results.

A last comment.  We have not tried to reference every relevant publication on the subjects of this work.  Such references are generally well known by active researchers in these areas and are readily accessible.  This view is consistent with the fact that the current article is not a review paper.

\section{Preliminaries}\label{Preliminaries}

Elementary heat conduction phenomena have been modeled by for approximately 150 years by the linear partial differential given in Eq.~(\ref{eq1}).  However, this PDE has associated with it two difficulties.  First,  Eq.~(\ref{eq1}) allows for the transformation of information between two separate space points with infinite speed (\cite{8}).  Second, this PDE does not provide a valid or accurate model for many heat conduction problems arising in a broad range of scientific and engineering situations (\cite{6,7,8}).  As a result of the existence of these issues,  there has been a large research effort to resolve and/or explain these issues; see for example (\cite{2, 5,6,7,8}).

In general,  the modifications to Eq.~(\ref{eq1}) have included adding a second-derivative term in time. The simplest such model is the Maxwell-Cattaneo representation (\cite{11, 4, 16}),
\begin{equation} \label{eq4}
\tau u_{tt} +u_t = Du_{xx}
\end{equation}
where the constant,  positive parameter $\tau$ is given the interpretation of a delay or lag time. Note that this PDE has exactly the same mathematical structure as the damped wave equation (\cite{2,7,11}).  Note that since the dependent variable,  $u(x,t)$,  represents the temperature as measured with respect to the absolute Kelvin scale, then only solutions satisfying the condition
\begin{equation} \label{eq5}
u(x,t) \geq 0, \quad t\geq 0,
\end{equation}
are physically meaningful. Inspection of Eq.~(\ref{eq4}) shows that it is a hyperbolic PDE and,  as a consequence,  can produce negative valued solutions.  So this is a difficulty that must be overcome.  However,  Eq.~(\ref{eq4}) does resolve the issue of infinite speed for information transfer.  We expect this speed to be approximately
\begin{equation} \label{eq6}
c=\sqrt{\frac{D}{\tau}}.
\end{equation}
Other modifications of Eq.~(\ref{eq1}) produce the linear PDE
\begin{equation*}
\tau u_{tt} +u_t=Du_{xx}+D_1u_{xxt},
\end{equation*}
where $(\tau,  D,  D_1)$ are non-negative,  constant parameters; see,  for examples, the paper of \cite{7}.  Here,  the issues also include whether the data can actually be measured to obtain the appropriate initial and boundary conditions.

Our analysis of these modified heat equations leads us to the conclusion that these efforts have no fundamental groundings in an overall consistent physical theory.  Further, they are related to other efforts involving the examination of macroscopic limits of microscopic mathematical models (\cite{1,3}) and the validity of expansions in a (supposedly) small parameter (\cite{12}).  An accurate,  but colorful way of characterizing this situation is to look at it from the perspective of a ``wack-a-mole" situation where the things required to resolve one issue provides the opportunity for the creation of another unresolved issue.

\section{Modeling and dynamic consistency}\label{Modeling}

A mathematical model of a system is a symbolic representation incorporating some aspects of its properties and features.  These aspects concern properties and important features of immediate interest to the modeler and the relevant scientific community.  However,  it must not be forgotten that a mathematical model,  in general,  can never capture all aspects and details of a system.  This is because at a given moment, no realistic system has all of its properties known and as a consequence,  the unknown ones can not  be included in the model.  Further,  even if a particular feature or aspect of a system is known,  the currently available mathematical structures may not be able to incorporate them in an appropriate and viable model.

In the next Section we present arguments which form the bases for the construction of a discrete microscopic model of heat conduction in one-space dimension.  Our dependent variable,  $u$,  represents the temperature measured on the absolute Kelvin scale.  This model will be derived by using the principle of dynamic consistency (\cite{13,15}).  We now provide a concise definition of dynamic consistency and explain how it may be applied to the construction of mathematical models.

\begin{definition}
Consider two systems, $S_1$ and $S_2$.  Let $S_1$ have some property, $P$.  If $S_2$ also has property $P$, then $S_2$ is said to be dynamically consistent to $S_1$ with respect to property $P$.
\end{definition}

Observe that the systems can be essentially anything. An explicit example is where $S_1$ is a differential equation and $S_2$ is a discretization of $S_1$.  The case where $S_2$ is a finite-difference discretization has been extensively studied within the context of the ``nonstandard finite differences” methodology (\cite{13,15}).

Another example is where $S_1$ is a physical system and $S_2$ is a mathematical model of $S_1$. For this case, $P$ might include any or all of the following features (\cite{13}):
\begin{enumerate}[(i)]
\item the dependent variables are non-negative
\item the dependent variables are bounded
\item conservation laws exist
\item the dependent variables have special asymptotic behaviors.
\end{enumerate}
Note that the mathematical model may not have all of the properties $P$. This is a general consequence of the modeling process. Again, considering the situation where $S_1$ is a differential equation and $S_2$ its discretization, then numerical instabilities may occur in the equations of discretizations (\cite{13,15}).

Since the major goals of mathematical modeling are to provide an understanding of the system, help to make valid predictions of the evolution of the physical system, and otherwise provide insights into the fundamental nature of the system,  we must be thoughtful in our construction of such models.  In particular,  we must always be aware that the results of mathematical modeling depend strongly on the mathematical structures we use to create the models.  An insightful example is the modeling of single-population systems where discrete models give solution behaviors entirely different from those provided by differential equations (\cite{14}).

Another difficulty in creating mathematical models is to take as exact particular functional forms which are based on experimental results valid only at low concentrations, or a limited range of temperatures, etc.  A prime illustration of this is Fick’s law for heat conduction (\cite{6,8}).

Finally,  we give a concise summary of the dynamic consistency methodology for constructing mathematical models of a system:
\begin{enumerate}[(a)]
\item Understand and analyze the system as fully possible;
\item from (a),  list the major properties and constraints that should be incorporated into the model;
\item using a particular mathematical structure (difference equations, differential equations, integral equations, etc.),  construct a mathematical model that includes as many as possible of the features and/or items listed in (b).
\end{enumerate}
Note that different mathematical models will arise from the effort to include some, all, or a separate set of properties of the system and its associated constraints,  such as symmetries and conservation laws.  In general,  the application of the dynamic consistency methodology does not yield a unique mathematical model (\cite{9,13}],  and often this ambiguity can be used to our advantage.  It should be clear that the inclusion of an increasing number of dynamic, consistent properties should yield better mathematical models.  However,  the modeler may not,  in general, know all the important properties of a system.  This implies that it is unlikely that the modeling process will generate equations that are ``exact,” either mathematically or physically.

In the next section,  we apply the dynamic consistency methodology to the issue of constructing a discrete microscopic level model for simple heat conduction.

\section{A discrete microscopic level model}\label{Discrete}

Before proceeding to construct the microscopic level model equation,  we need to provide some brief comments on several topics.

First,  we assume that it makes sense to divide the physical universe into three separate relevant domains in terms of space-time characterization.  This is depicted in Fig. ~(\ref{fig1}) and Table~(\ref{tbl1}).  The smallest or lowest level is where the atomic structure of matter dominates our understanding of physical systems. The top level,  sometimes called the continuum limit, is where macro-scale phenomena is often observed,  measured,  and interpreted.  Generally,  at this level,  the graininess of the atomic aspects of matter does not play an important role. Finally, the meso-level corresponds to considering an a priori atomic system,  but averaged over its atomic distinctness to produce variables that have macroscopic meaning.

Second,  we associate with each of the three levels characteristic time and space physical scales; see Fig. ~(\ref{fig1}). and Table~(\ref{tbl1}).  Thus, $(t_a, x_a)$ are constant parameters associated with the evolution of the system dynamics at the micro-level; $(\delta_t,\delta_x)$ are the constant parameters associated with physical phenomena at the meso-level; and $(T^*, L^*)$ are the time and length scales for what occurs at the macro-level.  Relations between their magnitudes are indicated in Fig.~(\ref{fig1}) and Table~(\ref{tbl1}).

\begin{figure}
\centering
\includegraphics[scale=0.1]{fig1.jpg}
\caption{Ordering of time and space scales}\label{fig1}
\end{figure}

\begin{table
\caption{Definition of time and space related variables}\label{tbl1}
\begin{tabular}{| l | l | l}
$u_s^r$  & $u_m^k$ & $u(x,t)$  \\
$t_r=t_a\cdot r$& $t_k = \delta_t\cdot k$ & $(x,t) \textrm{ continuous}$\\
$x_s = x_a\cdot s$ & $x_m = \delta_x \cdot m$ & \\
$r=0,1,2,\dots$ & $k=0,1,2, \dots$ & \\
$s= \textrm{ integers}$ & $m=\textrm{ integers}$ & \\
\end{tabular}
\end{table}

Our main purpose is to construct a discrete space-time model for the micro-level and then extend it to the meso-level.  Next,  we examine the meso-level model and interpret it as a modified or corrected version of what is usually considered to be a mathematical model for the continuum limit.

The space and time variables at the micro-level are discrete and are defined in Fig.~(\ref{fig1}) and Table~(\ref{tbl1}).  We show that the derived meso-level evolution equation has to be considered a hybrid differential-difference equation. This we explain below.  Note that $t_r$ and $x_s$ are discrete time and space variables.  The discrete meso-scale variables and their relationship to the continuum variables $(t,x)$ will also be clarified later.

To proceed in the derivative of a discrete micro-level model for heat conduction in one-dimension, we assume that our mathematical model satisfies the following conditions:
\begin{enumerate}[(A)]
\item be a linear difference equation,
\item satisfy the constraint that information can only travel with a finite speed in the discrete space-time;
\item be invariant under the parity transformation; i.e., the discrete equation should not change its form under the interchange
\begin{equation}\label{eq7}
(s+i) \ \longleftrightarrow \ (s-i),\\
s = \textrm{ integer}, \ i=1, 2, 3, \dots.
\end{equation}
\item not be invariant under time inversion;
\item satisfy a non-negativity condition for its solutions;
\item possess the property of measurability.\\
\end{enumerate}

The following is a brief discussion of the implications of these requirements:
\begin{enumerate}
\item[(A$_1$)] Assuming that the micro-level equation is linear provides a great simplification in that there are many more possible nonlinear equations and, generally, we have no clue as to which one(s) to select.
\item[(B$_1$)]  The micro-level is characterized by the time scale to $t_a$ and the space scale to $x_a$. Therefore, the speed of information flow,  $c$,  will be
\begin{equation} \label{eq8}
c = \frac{x_a}{t_a}
\end{equation}
\item[(C$_1$)]  For a discrete equation the condition given by Eq.~(\ref{eq7}) ensures parity conservation (\cite{13}).
\item[(D$_1$)]  The simplest way to have this condition hold is to have the difference equation depend only on the discrete times $t_r$ and $t_{r+1}$. Physically, this requirement implies physical dissipation.
\item[(E$_1$)]  Since our dependent variable corresponds to temperature on the absolution Kelvin scale, it must be always non-negative (\cite{6,11}).
\item[(F$_1$)]  By having the property of measurability, we mean that the model should not contain any parameter and/or dependent variables that can not be readily measured using standard experimental techniques.
\end{enumerate}

Let $u_s^r$ be the temperature at discrete time $t_r$ and location $x_s$.  The simplest difference equation that satisfies all of the above requirements is given by the expression
\begin{equation}
\begin{cases}\label{eq9}
u_s^{r+1} = &pu_{s+1}^r+(1-2p)u_s^r+pu_{s-1}^r\\
&0<p \leq \frac{1}{2}.
\end{cases}
\end{equation}
In the work to follow, we specialize to the case $p=1/3$ and Eq.~(\ref{eq9}) reduces to
\begin{equation} \label{eq10}
u_s^{r+1} = \left(\frac{1}{3}\right)\left(u_{s+1}^r+u_s^r+u_{s-1}^r\right).
\end{equation}
This special value of $p$ does not change any of our conclusion.

It should be indicated that Eq.~(\ref{eq10}) is a partial difference equation, first-order in the discrete-time and second-order in the discrete space variables.  Furthermore,  Eq.~(\ref{eq10}) can be solved exactly by a number of standard techniques (\cite{14}).  We will not provide such solutions since their explicit expressions play no role in the purpose of this article,  namely,  the derivation of the meso-level equation and how it should be interpreted.  However,  we do assume that there exists an analytic function, $U(x,t)$, such that
\begin{equation} \label{eq11}
\begin{cases}
u_s^r = U(x_s, t^r),\\
s=\textrm{ integer }; \quad r=(0,2,2,\dots); \quad x_a \textrm{ and } t_a \textrm{ are fixed.}
\end{cases}
\end{equation}
Note that Eq.~(\ref{eq10}) can be rewritten to the form
\begin{equation} \label{eq12}
u_s^{r+1}-u_s^r = \left(\frac{1}{3}\right)\left(u_{s+1}^r-2u_s^r+u_{s-1}^r\right),
\end{equation}
and from this we can obtain the expression
\begin{equation} \label{eq13}
\frac{u_s^{r+1}-u_s^r}{t_a} = D\left(\frac{u_{s+1}^r-2u_s^r+u_{s-1}^r}{x_a^2}\right),
\end{equation}
where
\begin{equation} \label{eq14}
D=\left(\frac{1}{3}\right)\left(\frac{x_a^2}{t_a}\right).
\end{equation}
If $D$ is assumed constant and if the limits
\begin{equation}\label{eq15}
t_a \rightarrow 0,  \quad x_a \rightarrow 0,  \quad D = \textrm{ constant},
\end{equation}
then Eq.~(\ref{eq13}) becomes the following partial differential equation
\begin{equation} \label{eq16}
u_t=Du_{xx}.
\end{equation}

\begin{comm}
In more detail, the following are the correct limits:
\begin{equation}
\begin{cases}
t_a \rightarrow 0, x_a \rightarrow 0; \quad  r \to \infty,  \quad s \to \infty;\\
\left(\frac{1}{3}\right)\left(\frac{x_a^2}{t_a}\right) = D = \textrm{ constant};\\
t_ar = t = \textrm{ fixed},  \quad x_as = x = \textrm{ fixed.}
\end{cases}
\end{equation}
\end{comm}

Observe that Eq.~(\ref{eq16}) is the standard heat equation (\cite{8}).  Also, the procedure used to derive this result opens up a number of important issues related to the macroscopic limits of microscopic models (\cite{1,3,9}).  As we shall see in the next section,  an alternative analysis and interpretation is possible.

\section{Macro-level approximation}

Consider Eq.~(\ref{eq10}) and as stated above,  assume that there exists a function, $U(x,t)$,  with appropriate analytical properties in $x$ and $t$,  such that
\begin{equation} \label{eq18}
u_s^r=U(x_x,t_r).
\end{equation}
This implies that $U(x,t)$, for fixed $\overline{x}$ and $\overline{t}$,  satisfies the same recurrence/different equation as $u_s^r$, i.e.,
\begin{equation} \label{eq19}
U(x,t+\overline{t})=\left(\frac{1}{3}\right)\left[U(x+\overline{x},t)+U(x,t)+U(x-\overline{x},t)\right]
\end{equation}
Thus, at the meso-level where $(\delta_x,\delta_t)$ are the length and time (physical) scales,  we may identify $(\overline{x},\overline{t})$, respectively, with $(\delta_x,\delta_t)$,  and obtain
\begin{equation}\label{eq20}
U(x,t+\delta_t)=\left(\frac{1}{3}\right)\left[U(x+\delta_x,t)+U(x,t)+U(x-\delta_x,t)\right]
\end{equation}
To proceed,  carry out Taylor series expansions of the individual terms in Eq.~(\ref{eq20}) and retain only terms up to $O(\delta_t^3)+O(\delta_x^6)$.  Doing this gives, after division by $\delta_t$, the expression
\begin{equation}\label{eq21}
U_t+\left(\frac{\delta_t}{2}\right)U_{tt}+O(\delta_t^2) = DU_{xx}+\left[\left(\frac{D}{36}\right)\delta_x^2\right]U_{xxxx}+O(D \delta_x^4),
\end{equation}
where
\begin{equation}\label{eq22}
D=\frac{\delta_x^2}{3\delta_t}.
\end{equation}
Note by assuming that $D$ is fixed in value,  we have
\begin{equation}\label{eq23}
\delta_t=O(\delta_x^2).
\end{equation}
Consequently,  to terms of $O(1)$ and $O(\delta_t)+O(\delta_x^2)$, the evolution equations for $U(x,t)$ are
\begin{align}
O(1) &: \overline{U}_t=D\overline{U}_{xx},\label{eq24} \\
O(1)+\mathcal{O}(\delta_t)+O(\delta_x^2) &: \left(\frac{\delta_t}{2}\right)\overline{U}_{tt}+\overline{U}_t=D\overline{U}_{xx}+\left[\left(\frac{D}{36}\right)\delta_x^2\right]\overline{U}_{xxxx}. \label{eq25}
\end{align}
We have placed bars over the $U$’s to indicate that these equations only yield approximations to the exact PDE for $U(x,t)$.

Note that the lowest order expansion,  i.e.,  retaining terms to $O(1)$, the evolution equation is the standard heat conduction PDE.  Since
\begin{equation}\label{eq26}
\partial_t=D\partial_{xx}+O(\delta_t)+O(\delta_x^2),
\end{equation}
we have
\begin{align}\label{eq27}
\partial_{xxxx} = &\partial_{xx}\cdot \partial_{xx}\cr
= &\partial_{xx}\left[\left(\frac{1}{D}\right)\partial_t\right]+O(\delta_t)+O(\delta_x^2).
\end{align}
Substituting this latter result into the second term on the right side of Eq.~(\ref{eq25}) gives
\begin{equation} \label{eq28}
\left(\frac{\delta_t}{2}\right)\overline{U}_{tt}+\overline{U}_t=D\overline{U}_xx+\left(\frac{\delta_x^2D}{36}\right)\overline{U}_{xxt}.
\end{equation}
This PDE is the same PDE that was derived by \cite{7} under different assumptions.

\section{Discussion}

In this work, we began by assuming that we have the exact micro-level evolution equation for simple heat conduction,  given by Eq.~(\ref{eq10}).  Assuming that the solution of this discrete equation can be represented as the values of an analytic function,  $U(x,t)$,  evaluated on the lattice of discrete space and time points,  we used Taylor series expansion on this expression at the meso-level to obtain what we will call modified hybrid-PDEs,  see for example Eq.~(\ref{eq28}).  Dropping the bar over $\overline{U}$,  the following two lowest order PDEs were found
\begin{equation} \label{eq29}
U_t=DU_{xx},
\end{equation}
\begin{subequations}
\vspace{-.3in}
\begin{align}
\tau U_{tt}+U_t = &DU_{xx}+D_1U_{xxt},  \label{eq30a}\\
\tau U_{tt}+U_t = &DU_{xx}+D_2U_{xxxx},  \label{eq30b}
\end{align}
\end{subequations}
where in this representation, we have expressed them, especially the Eqs.~(\ref{eq30a}, \ref{eq30b}),  in a general way suitable for use as phenomenological equations, where $(\tau,D,D_1)$ are to be determined from experiment. In principle, as indicated in the previous section, these parameters are determined by the time and length scales at the meso-level.

Observe that Eq.~(\ref{eq29}) is the standard heat equation. Within our interpretation, this PDE should be considered the macroscopic limiting equation of the microscopic Eq.~(\ref{eq10}).  We do not designate it by the name continuum limit (\cite{1,3}), since no micro-level or meso-level parameters are required to go to zero.  Also, note that within the context of this work, relationships should exist between the space and time scales at the micro- and meso-levels, e.g.,
\begin{equation} \label{eq31}
\delta_x=N_1x_a, \quad \delta_t=N_2t_a,
\end{equation}
where $N_1$ and $N_2$ are large integers.

Since neither the micro- or meso-levels parameters go to zero, we expect the small parameters $(\epsilon_1, \epsilon_2)$, where
\begin{equation}\label{eq32}
\epsilon_2=\frac{\delta_t}{T^*}, \quad \epsilon = \frac{\delta_x}{L^*},
\end{equation}
to play important roles in the expression of the first modification to the macroscopic equation.  To show this, carry out the scaling of independent variables in Eq.~(\ref{eq28}) as follows:
\begin{equation}\label{eq33}
t= T^*\overline{t}, \quad x=L^*\overline{x},
\end{equation}
where $(L^*, T^*)$ are the macro-level scales,  and $(\overline{x}, \overline{t})$ are the dimensionless independent variables.  Since the PDE is linear,  we need to make no explicit scaling for $U$. Substituting Eq.~(\ref{eq31}) into Eq.~(\ref{eq28}) and multiplying each term by $T^*$ gives
\begin{equation} \label{eq34}
\left[ \left( \frac{1}{2} \right) \left( \frac{\delta_t}{T^*}\right) \right] \frac{\partial^2 U}{\partial \bar{t}^2}+\frac{\partial U}{\partial \bar{t}}
=\left( \frac{DT^*}{(L^*)^2}\right)\frac{\partial^2U}{\partial \bar{x}^2}+\left( \frac{\delta_x^2}{36}\right)\left[ \frac{T^*D}{(L^*)^4} \right]\ \frac{\partial^4U}{\partial \bar{x}^4}.
\end{equation}
Keeping in mind the definition of $D$ given in Eq.~(\ref{eq22}), we define $\overline{D}$,  $\epsilon_1$, and $\epsilon_2$ as follows
\begin{equation} \label{eq35}
\overline{D}=\frac{DT^*}{(L^*)^2}, \quad \epsilon_1=\frac{1}{2}\left( \frac{\delta_t}{T^*} \right),  \quad \epsilon_2=\left( \frac{\delta_x}{L^*}\right)^2,
\end{equation}
where it should be noted that these three parameters are dimensionless,  i.e.,  pure numbers.  Also, it is expected that $\epsilon_1$ and $\epsilon_2$ are small,  i.e.,
\begin{equation}\label{eq36}
0<\epsilon_1<<1, \quad 0 < \epsilon_2<<1.
\end{equation}
With this notation,  Eq.~(\ref{eq28}) in dimensionless form is
\begin{equation} \label{eq37}
\epsilon_1 U_{\bar{t} \,\bar{t}}+U_{\bar{t}}=\overline{D}U_{\bar{x}\bar{x}}+\epsilon_2\overline{D}U_{\bar{x}\bar{x}\bar{x}\bar{x}}.
\end{equation}
From the fact that $\delta_t=O(\delta_x^2)$, we can expect an infinite hierarchy of modified PDEs,  starting from
\begin{equation} \label{eq38}
\begin{cases}
O(1): U_t=DU_{xx},\\
O(\delta_t)+O(\delta_x^2):  \textrm{Eq.~(\ref{eq37})}
\end{cases}
\end{equation}

\begin{figure
\centering
\includegraphics[scale=0.1]{Fig2a.jpg}
\caption{Hiearchy of modified PDE.}\label{fig2}
\end{figure}

Keep in mind the fact the $N$-th modified PDE contains terms from all the previous equations. The $N=0$ and $N=1$ cases are the two expressions presented in Eq.~(\ref{eq38}).

Observe that at the $N$-th level of the modified PDEs hierarchy,  the time derivatives up to order $(N+1)$ appear.  This means also that spatial derivatives to even order $2(N+1)$ must occur.  Thus,  for $N=1$,  the time derivatives $U_t$ and $U_{tt}$ appear only along with the spatial derivatives, $U_{xx}$ and $U_{xxxx}$ or $U_{xxt})$.  However,  for most if not all one-dimension space systems,  the initial condition,  $U(x,0)$,  is known,  along with related boundary conditions, but $U_t(x,0)$ or higher time derivative are not known, it follows that none of the modified PDEs is measurable,  i.e.,  their mathematical solutions require information that can not be physically determined.  An important implication is that this way or methodology of investigating heat condition is invalid or at least not productive.  One way to view this situation is to remember that if we begin with a discrete space and time model, then the mathematical structure currently used corresponds to difference equations.  However,  it is well known that difference equations do not model exactly differential equations; see the books by (\cite{14,15}).  In particular,  the following relation holds
\begin{equation} \label{eq39}
\textrm{ difference equation } = \textrm{ infinite-order differential equation.}
\end{equation}
This can easily be seen by considering the following argument:
\begin{itemize}
\item[(1)] Define the (difference) shift and derivative operators by the following relations
\begin{equation}\label{eq40}
E_af(x) \equiv f(x+a), \quad Df(x)\equiv \frac{df(x)}{dx}.
\end{equation}
\item[(2)] The Taylor expansion of $f(x+a)$ is
\begin{align}\label{eq41}
f(x+a)&=\overset{\infty}{\underset{k=1}{\sum}} \left(\frac{a^k}{k!}\right)D^kf(x)\cr
&=e^{aD}f(x).
\end{align}
\item[(3)] From the definition of the shift operator,  it follows that the operators satisfy
\begin{equation} \label{eq42}
E_a=e^{aD}
\end{equation}
\item[(4)] The difference operator is defined by the relation
\begin{equation} \label{eq43}
\Delta_a=E_a-1.
\end{equation}
Comparison of Eq.~(\ref{eq42}) and Eq.~(\ref{eq43}) gives
\begin{equation} \label{eq44}
\Delta_a=e^{aD}-1.
\end{equation}
\item[(5)] Since $\Delta_a/a$ and $D$,  are the discrete and continuous forms of the derivative, Eq.~(\ref{eq44}) shows that they nonlinearly related and any finite order difference equation corresponds to an infinite order differential equation.  Also,  if we solve Eq.~(\ref{eq44}) for $D$, we obtain
\begin{align}
 D &\equiv \left(\frac{1}{a}\right) \ln (1+\Delta_a)\\
&= \left(\frac{1}{a}\right)\left[\Delta_a-\frac{\Delta_a^2}{2}+\frac{\Delta_a^3}{3}-\frac{\Delta_a^4}{4}+\cdots\right].
\end{align}
Consequently,  it can be concluded that in general a finite order differential equation is equivalent to an infinite order difference equation.
\end{itemize}

In summary,  we have used the dynamic consistency methodology (\cite{13,14}) to construct a discrete space and time microscopic mathematical model for one-space dimension heat conduction.  To obtain solutions to this difference equation, only initial and/or boundary values are required to be known. After some manipulations,  we were able to derive an infinite number of mesoscopic PDEs.  These PDEs form a hierarchy of equations,  each of which modifies the equations on previous levels; see Fig. ~(\ref{fig2})  However, the modified equations at the first and higher levels generally require knowledge of time derivatives for which experimental data is either difficult or impossible to measure.  Our major conclusion is that the procedures used to construct and justify previous generalizations of the basic heat conduction equation,  $u_t=Du_{xx}$,  are not valid and that other mathematical structures must be discovered or created to resolve issues such as infinite speed of propagation for the transfer of information (\cite{2,4,5,6,7,8,11,16}).  We are currently investigating this problem.

\section{Bibliography}\label{Bibliography}

\bibliographystyle{cas-model2-names}


\section*{Abstract}
{\small

After a long qualifying process packed with surprises (Italy missing out as the reigning European champions) and last minute drama (both Egypt and Peru missed out on penalties), the FIFA World Cup 2022 kicked off on the 20th of November in Qatar. With 32 countries and over 800 players representing nearly 300 clubs globally, it measured up to more than 12 billion EUR in the players’ current estimated market value total. In this short piece, we explore what the small and interconnected world of football stars looks like and even make a few efforts and compare success in soccer to social networks.

}

\vspace{0.5cm}
{\small {\bf Keywords}: network science, social network analysis, soccer, data science}

\vspace{1.0cm}
{\it \hspace{-1cm} Published in Nightingale, Journal of the Data Visualization Society, December 23, 2022~\cite{nightingale}. }
\vspace{1.0cm}

\section{Data}

We are data scientists with a seasoned football expert on board, so we went for one of the most obvious choices of the field – www.transfermarkt.com. We first wrote a few lines of Python code to scrape the list of participating teams~\cite{teams}, the list of each team’s players~\cite{players}, and the detailed club-level transfer histories of these players arriving at the impressive stats of our intro by comprising the complete transfer history of 800 players, measuring up to 6,600 transfers and dating back to 1995 with the first events.

\section{Club network}

The majority of players came from the top five leagues (England, Spain, Italy, Germany, and France) and represented household teams such as Barcelona (with 17 players), Bayern Munich (16), or Manchester City (16). While that was no surprise, one of the many wonders of a World Cup is that players from all around the globe can show their talents. Though not as famous as the ‘big clubs’, Qatari Al Sadd gave 15 players, more than the likes of Real Madrid or Paris Saint-Germain! There are, however, great imbalances when throwing these players’ market values and transfer fees into the mix. To outline these, we decided to visualize the typical ‘migration’ path football players follow – what are the most likely career steps they make one after the other?

A good way to capture this, following the prestige analysis of art institutions~\cite{fraiberger2018quantifying}, is to introduce network science~\cite{netsci} and build a network of football clubs. In this network, every node corresponds to a club, while the network connections encode various relationships between them. These relationships may encode the interplay of different properties of clubs, where looking at the exchange of players (and cash) seems a natural choice. In other words, the directed transfers of players between clubs tie the clubs into a hidden network. Due to its directness, this network also encodes information about the typical pathways of players via the ‘from’ and ‘to’ directions, which eventually capture the different roles of clubs as attractors and sinks.

\begin{table}[!hbt]
\centering
  \includegraphics[scale=0.5]{table1.png}
  \caption{The datafied transfer history of Neymar.}
  \label{tab:tab1}
\end{table}

To do this in practice, our unit of measure is the individual transfer history of each player, shown in Table \ref{tab:tab1} for the famous Brazilian player known simply as Neymar. This table visualizes his career trajectory in a datafied format, attaching dates and market values to each occasion he changed teams. His career path looks clean from a data perspective, although football fans will remember that it was anything but – his fee of EUR 222M from Barcelona to PSG still holds the transfer record to this day. These career steps, quantified by the transfers, encode upgrades in the case of Neymar. In less fortunate situations, these prices can go down signaling a downgrade in a player’s career.

\begin{figure}[!hbt]
\centering
\includegraphics[width=0.75\textwidth]{Figure1.png}
\caption{The network of the top football clubs based on the total amount of money spent and received on player transfers. Node sizes correspond to these values, while node coloring shows the dominant color of each club’s home country flag.}
\label{fig:fig1}
\end{figure}

\clearpage

Following this logic in our analysis, we assumed that two clubs, A and B, were linked (the old and new teams of a player), if a player was transferred between them, and the strength of this link corresponded to the total amount of cash associated with that transaction. The more transactions the two clubs had, the stronger their direct connection was (which can go both ways), with a weight equal to the total sum of transfers (in each direction). In the case of Neymar, this definition resulted in a direct network link pointing from Barcelona to Paris SG with a total value of EUR 222M paid for the left winger.

Next, we processed the more than six thousand transfers of the 800+ players and arrived at the network of teams shown in Figure \ref{fig:fig1}. To design the final network, we went for the core of big money transactions and only kept network links that represented transfer deals worth more than EUR 2.5M in total. This network shows about 80 clubs and 160 migration channels of transfers. To accurately represent the two aspects of transfers (spending and earning) we created two versions of the same network. The first version measures node sizes as the total money invested in new players (dubbed as spenders), while the second version scales nodes as the total money acquired by selling players (dubbed as mentors).

\paragraph{Spenders.}

The first network shows us which clubs spent the most on players competing in the World Cup, with the node sizes corresponding to the total money spent. You can see the usual suspects: PSG, the two clubs from Manchester, United, and City, and the Spanish giants, Barcelona, and Real Madrid. Following closely behind are Chelsea, Juventus, and Liverpool. It’s interesting to see Arsenal, who – under Arteta’s management – can finally spend on players, and Bayern Munich, who spend a lot of money but also make sure to snatch up free agents as much as possible.

Explore these relationships and the network in more detail by looking at Real Madrid! Los Blancos, as they’re called, have multiple strong connections. Their relationship with Tottenham is entirely down to two players who played an integral part in Real Madrid’s incredible 3-year winning spell in the Champions League between 2016 and 2018: Croatian Luka Modric cost 35M, and Welsh Gareth Bale cost an at-the-time record-breaking 101M. While Real Madrid paid 94M for Cristiano Ronaldo in 2009 to Man Utd, in recent years there was a turn in money flow, and United paid a combined 186M for three players: Ángel Di María, Raphael Varane, and Casemiro. They also managed to sell Cristiano Ronaldo with a profit to Juventus for 117M.

One can see other strong connections as well, such as Paris SG paying a fortune to Barcelona for Neymar and Monaco for Kylian Mbappé. There are also a few typical paths players take – Borussia Dortmund to Bayern Munich, Atlético Madrid to Barcelona, or vice versa. It’s also interesting to see how many different edges connect to these giants. Man City has been doing business worth over EUR 1M with 27 different clubs.

\paragraph{Mentors.}

The second network shows which clubs grow talent instead of buying them and have received a substantial amount of money in return. Node sizes represent the amount of transfer fees received. This paints a very different picture from our first network except for one huge similarity: Real Madrid. In the past, they were considered the biggest spenders. They have since adopted a more business-focused strategy and managed to sell players for high fees as mentioned above.

A striking difference, however, is while the top spenders were all part of the top five leagues, the largest talent pools came from outside this cohort except for Monaco. Benfica, Sporting, and FC Porto from Portugal, and Ajax from the Netherlands are all famous for their young home-grown talents, and used as a stepping stone for players from other continents. Ajax has sold players who competed in this World Cup for over EUR 560M. Their highest received transfer fees include 85.5M for Matthijs de Ligt from Juventus and 86M for Frenkie de Jong from Barcelona. Ajax signed de Jong for a total of EUR 1 from Willem II in 2015 when he was 18, and de Light grew up in Ajax’s famous academy. Not to mention that they recently sold Brazilian Antony to Manchester United for a record fee of 95M. They paid 15.75M for him just 2 years ago – that’s almost 80M in profit. Insane!

Benfica earned close to 500M, most recently selling Uruguayan Darwin Nunez for 80M to Liverpool. The record fee they received is a staggering 127M for Portuguese Joao Félix from Atlético Madrid, who grew up at Benfica. Monaco earned 440M from selling players such as Kylian Mbappé (180M) and Aurélien Tchouameni (80M), Portuguese Bernardo Silva (50M), Brazilian Fabinho and Belgian Youri Tielemans (both for 45M). These clubs have become incredible talent pools for the bigger clubs, therefore really appealing to young players. It’s interesting to see how many edges the nodes for these clubs have, further proving that these teams function as a means for reaching that next level.

\section{Player network}

After looking at the club-to-club relationships, zoom in on the network of players binding these top clubs together. Here, we built on the players’ transfer histories again and reconstructed their career timelines. Then we compared these timelines between each pair of World Cup players, noted if they ever played for the same team, and if so, how many years of overlap they had (if any).

To our biggest surprise, we got a rather intertwined network of 830 players connected by about 6,400 former and current teammate relationships, as shown in Figure \ref{fig:fig2}. Additionally, the so-called average path length turned out to be 3 – which means if we pick two players at random, they most likely both have teammates who played together at some point. Node sizes were determined by a player’s current market value, and clusters were colored by the league’s nation where these players play. It didn’t come as a surprise that current teammates would be closer to each other in our network. You can see some interesting clusters here, with Real Madrid, Barcelona, PSG, and Bayern Munich dominating the lower part of the network and making up its center of gravity.

Why is that? The most valuable player of the World Cup was Kylian Mbappé, with a market value of 160M, surrounded by his PSG teammates like Brazilians Marquinhos and Neymar and Argentinian Lionel Messi. Messi played in Barcelona until 2021, with both Neymar and Ousmane Dembélé connecting the two clusters strongly. Kingsley Coman joined Bayern Munich in 2017, but he played for PSG up until 2014, where they were teammates with Marquinhos, thus connecting the two clusters.

You can discover more interesting patterns in this graph, such as how the majority of the most valuable players have played together directly or indirectly. You can also see Englishmen Trent Alexander-Arnold (Liverpool) or Declan Rice (West Ham United) further away from the others. Both of those players only ever played for their childhood clubs. But the tight interconnectedness of this network is also evident with how close Alexander-Arnold actually is to Kylian Mbappé. During the 2017–2018 season, Mbappé played at Monaco with Fabinho behind him in midfield, who signed for Liverpool at the end of the season, making him and Alexander-Arnold teammates.

\clearpage
\begin{figure}[!hbt]
\centering
\includegraphics[width=0.90\textwidth]{Figure2.png}
\caption{The player-level network showing previous and current teammate relationships. Note size corresponds to the players’ current market values, while color encodes their nationality based on their country’s flag’s primary color. See the interactive version of this network here~\cite{interactive}.}
\label{fig:fig2}
\end{figure}

With the World Cup hosting hundreds of teams’ players from various nations, there are obviously some clusters that won’t connect to these bigger groups. Many nations have players who have only played in their home league, such as this World Cup’s host nation Qatar (maroon cluster in the top left corner). Saudi Arabia (green cluster next to Qatar) beat Argentina, causing one of this year’s biggest surprises. Morocco (red cluster in the top right corner) delivered the best-ever performance by an African nation in the history of the World Cups. Both of those nations join Qatar in this category of home-grown talent. These players will only show connections if they play in the same team – in the case of the Moroccan cluster, that team is Wydad Casablanca. The Hungarian first league’s only representative at the World Cup, Tunisian Aissa Laidouni from Ferencváros hasn’t played with anyone else on a club level who has made it to the World Cup. He became a lone node on our network. That shouldn’t be the case for long, considering how well he played in the group stages.

\section{Success and networks}

The potential success of different teams and the outcomes of championships have been majorly interesting for data and statistics people since the era of Moneyball. Ever since, a wide range of efforts came to light about the possibilities to predict the outcomes of sporting events, from asking actual animals like Paul the Octopus to serious academic research, and even companies specializing in this domain.

Here we are not attempting to overcome such elaborate methods and solutions but intend to show the quantitative drivers of success in soccer from a different angle - as network science sees it. For that, we follow some of the earlier work of Barabasi et. al and others on the quantification of success and the role of networks in team success~\cite{guimera2005team, barabasi2018formula, janosov2020quantifying}.

\paragraph{Player-level success.} To capture the success of players, one could pick a large number of different KPIs depending on the goal they wish to study. Here, to illustrate the role of networks, we picked their current market value as a proxy of their success. We also note that we do not take sides on how much the player's financial success may proxy their talent or their sense of branding, just expect what the market shows.

To assess the players' characteristics, we computed several measures describing them. First, we added measures describing their career: the year of their first transfer, their first public market value, and the number of transfers they had. Then we added several network measures capturing their network positions: Weighted Degree measuring the total weight of connections a player has, Closeness Centrality measuring how few hops a player is from the others, Betweenness Centrality centrality capturing the network bridge positions, and Clustering tells us what fraction of a player's connections know each other.

\begin{figure}[!hbt]
\centering
\includegraphics[width=0.80\textwidth]{Figure3.png}
\caption{The number of players, and the correlations between their market value and the different descriptive features as a function of the number of transfers they had.}
\label{fig:fig3}
\end{figure}

Then, we conducted a series of correlation analyses where we correlated these measures against each player's current market value. Additionally, we differentiated between players based on their seniority, measured by the number of transfers they had. We realize having played for more teams doesn't necessarily mean that a player has had a longer career, younger players can be loaned out to gain experience, but this approach works well with our data set. Figure \ref{fig:fig3} shows the effect of this on the number of players analyzed: less than 1\% of the players have less than 2 transfers, while the line between 7 and 8 transfers splits
the players into two roughly equal parts. Also, slightly less than 5\% of the players have more than 15 transfers which translates into about 30 veterans.

Furthermore, Figure \ref{fig:fig3} also shows the correlations between the players' descriptors against their seniority. We sorted this chart based on the full population of players, meaning that the features in the correlation matrix are put in a descending order based on the correlation matrix when the number of transfers equals zero. The changing patterns in the matrix indicate the different correlation trends as we restrict our analysis to more and more senior players.

The most striking feature of the correlation values, ranging from 0.53 to -0.62, is that no matter the player's seniority, Weighted Degree plays a superior role. One can also notice that the negative correlation values typically occur for the most senior players, and are probably more accounted for noise than signal. For younger players, the graph also shows that after the Weighted Degree comes to their first market value, and then Closeness Centrality. However, at the tenure of about five transfers, the picture changes, and the top three most correlated features to the current market value end up being derived from the network of players. While we are not claiming pinpoint accuracy here, these findings certainly indicate that the role of networks shows significant potential for the expected success of a soccer player.

\paragraph{Country-level success.} Last but not least, we also shoot our shot at understanding some aspects of the final - and certainly not unsurprising - ranking of the World Cup. For this, we used the final ranking of all 32 teams at the 2022 FIFA World Cup~\cite{sportingnews}. Then, we went back to the previously introduced player-level statistics covering the total number of transfers, the first year of transfer, the total market value, and several network measures first attached to each player, and then aggregated to the level of the country's teams. The aggregation, depending on the distribution of the underlying values, happened either by taking the mean or the median values of a country's team members.

Next, we ventured into the world of predicting algorithms. To be more precise, we built a binary XGBoost~\cite{chen2016xgboost} classifier aimed to distinguish between those teams which made it to the Rd of 16 versus those that didn't make it longer than the Groups. We note that even in this scenario, the number of data points was very low. In future research, this part of the analysis could be improved by adding other World Cups to the data, which may also allow one to do more elaborate predictions, covering the Quarterfinals or even the grand finale.

\begin{figure}[!hbt]
\centering
\includegraphics[width=0.7\textwidth]{Figure4.png}
\caption{Feature importance expressed by relative values when differentiating the teams tho dropped out after the first round from the rest by using binary classification.}
\label{fig:fig4}
\end{figure}

We optimized our XGBoost model using 250 estimators, a maximum tree depth of 6, a learning rate varying between 0.1 and 0.005, a grid search algorithm, and 5-fold cross-validation. On the one hand, the overall prediction accuracy resulted in a fairly modest value of 60\%. Still, it is worth taking a look at the feature importance analysis shown in Figure \ref{fig:fig4}. This figure shows that - probably not so - surprising, the most important driver of a team's position is the total current market value of its players. High paychecks seem to pay off. At least to some extent pays - the correlation between team rank and the total value is still just somewhere between 0.5 and 0.6. This is then followed, somewhat in a tie, between closeness centrality and first market value. While first market value - for more junior players - is certainly correlated with their current value, a network centrality reaching that level certainly emphasizes the role of the network. While the other network measures are performing around 10\% each in terms of their relative importance, we leave it for previous research to nail down the exact meaning of these figures.

Finally, to outline the signal behind the noise caused by the low number of data points, we created Figure \ref{fig:fig5} showing how closeness centrality and market value guide the teams in terms of their final rankings.

\begin{figure}[!hbt]
\centering
\includegraphics[width=0.5\textwidth]{Figure5.png}
\caption{On this scatter plot, each team is represented by a dot colored according to its final rank, while its position is determined by the median closeness centrality and the total current market value of its players.}
\label{fig:fig5}
\end{figure}

\section{Conclusion}

In conclusion, we saw in our analysis how network science and visualization can uncover and quantify things that experts may have a gut feeling about but lack the hard data. This depth of understanding of internal and team dynamics that is possible through network science can also be critical in designing successful and stable teams and partnerships. Moreover, this understanding can lead to exact applicable insights on transfer and drafting strategies or even spotting and predicting top talent at an early stage. While this example is about soccer, you could very much adapt these methods and principles to other collaborative domains that require complex teamwork and problem solving with well-defined goals, from creative production to IT product management.