# AdaSub: Stochastic Optimization Using Second-Order Information in Low-
Dimensional Subspaces††thanks: This work was supported by the Novo Nordisk
Foundation under grant number NNF20OC0061894.

João Victor Galvão da Mata  Department of Applied Mathematics and Computer
Science
Technical University of Denmark
<EMAIL_ADDRESS>
ORCID 0009-0009-8251-6476  Martin S. Andersen  Department of Applied
Mathematics and Computer Science
Technical University of Denmark
<EMAIL_ADDRESS>
ORCID 0000-0002-4654-3946

###### Abstract

We introduce AdaSub, a stochastic optimization algorithm that computes a
search direction based on second-order information in a low-dimensional
subspace that is defined adaptively based on available current and past
information. Compared to first-order methods, second-order methods exhibit
better convergence characteristics, but the need to compute the Hessian matrix
at each iteration results in excessive computational expenses, making them
impractical. To address this issue, our approach enables the management of
computational expenses and algorithm efficiency by enabling the selection of
the subspace dimension for the search. Our code is freely available on GitHub,
and our preliminary numerical results demonstrate that AdaSub surpasses
popular stochastic optimizers in terms of time and number of iterations
required to reach a given accuracy.

###### Index Terms:

Stochastic optimization, subspace optimization, low-dimensional optimization,
stochastic quasi-Newton methods.

##  I Introduction

The optimization of neural network (NN) models is often a challenging task due
to the high-dimensional and non-convex nature of many such problems. First-
order optimization methods such as stochastic gradient descent (SGD) and its
many variants are widely used as the main optimization methods for training NN
models [4, 5, 11], but they often require extensive hyperparameter tuning to
reach an acceptable training loss [16] and only rely on noisy first-order
information, resulting in slow convergence and poor generalization for ill-
conditioned training problems [1].

Hyperparameters, including learning rate, weight decay, and momentum
parameters, are essential for the training process and require careful tuning.
However, tuning these parameters can be costly in terms of both computational
resources and human effort. In addition, a poor choice of optimizer and/or
parameters can lead to significant performance degradation, highlighting the
limitations of first-order methods [1, 16].

Optimization methods that incorporate second-order information can more
effectively capture and leverage the curvature of the loss function. As a
result, they have the potential to find better search directions. However,
this improved performance generally comes at the cost of additional
computational resources.

This paper presents “AdaSub,” a stochastic optimizer that utilizes second-
order information in a low-dimensional subspace without requiring the explicit
computation of the Hessian matrix. As a result, our method avoids significant
increases in the computational cost. Our approach is inspired by recent
methods such as AdaHessian [16] and SketchySGD [6], which access second-order
information through Hessian-vector products. Specifically, given a twice
continuously differentiable function $f\colon\mathbb{R}^{n}\to\mathbb{R}$ and
vectors $v\in\mathbb{R}^{n}$ and $w\in\mathbb{R}^{n}$, we have that

$\left.\frac{d\nabla f(w+tv)}{dt}\right|_{t=0}=\frac{\partial(\nabla
f(w)^{T}v)}{\partial w}=\nabla^{2}f(w)v,$ (1)

which can be computed without forming $\nabla^{2}f(w)$ by leveraing automatic
differentiation. SketchySGD uses such Hessian-vector products to compute a
diagonal plus low-rank approximation to $\nabla^{2}f(w)$ using sketching
techniques, and AdaHessian uses similar ideas to compute a diagonal
approximation to $\nabla^{2}f(w)$. In contrast, we use Hessian-vector products
to compute a second-order approximation of the objective function $f$ on a
subspace that is defined adaptively based on a current (stochastic) gradient
and a small set of past (stochastic) gradients rather than random directions.
This paper presents the following contributions:

  1. 1.

We propose a new stochastic optimization algorithm that generalizes SGD by
incorporating second-order information in a low-dimensional subspace. The
dimension of the subspace is a parameter of the algorithm.

  2. 2.

We provide a basic convergence analysis of our method in the special case
where the objective function is strongly convex and smooth.

  3. 3.

We demonstrate the effectiveness of AdaSub through numerical experiments based
on various learning tasks. Our method matches or outperforms state-of-the-art
methods in terms of the number of iterations.

  4. 4.

Our preliminary implementation is available on GitHub111Available at:
https://github.com/Jvictormata/adasub.

##  II Adaptive Subspace Minimization

The basic idea behind our method is to first compute a second-order
approximation of a twice continuously differentiable objective function
$f\colon\mathbb{R}^{n}\to\mathbb{R}$ on an $m$-dimensional affine set. Given a
vector $w\in\mathbb{R}^{n}$ and a matrix $V\in\mathbb{R}^{n\times m}$ with
orthonormal columns and with $m\ll n$, the restriction of $f$ to the affine
set $w+\operatorname{Range}(V)$ may be parameterized as $h(z)=f(w+Vz)$ where
$z\in\mathbb{R}^{m}$. Using the chain rule, we can express the second-order
Taylor expansion of $h$ around $0$ as

$\displaystyle\tilde{h}(z)$ $\displaystyle=h(0)+\nabla
h(0)^{T}z+\frac{1}{2}z^{T}\nabla^{2}h(0)z$
$\displaystyle=f(w)+g^{T}Vz+\frac{1}{2}z^{T}V^{T}HVz,$ (2)

where we use the shorthand notation $g=\nabla f(w)$ and $H=\nabla^{2}f(w)$. In
the special case where $\nabla^{2}h(z)\succ 0$, the unique minimizer of
$\tilde{h}$ is given by $z^{\star}=-(V^{T}HV)^{-1}V^{T}g$, which motivates an
update of the form

$w^{+}\leftarrow w+\eta\,Vz^{\star}=w-\eta\,V(V^{T}HV)^{-1}V^{T}g,$ (3)

where $\eta>0$ is a step size. This corresponds to a Newton step for the
restriction of $f$. The computation of $z^{\star}$ requires $m$ Hessian-vector
products to form $HV$, and hence the cost of computing $z^{\star}$ grows as
$O(mC+m^{3})$ if $C$ is the cost of a Hessian-vector product. This is much
lower than the cost of computing a Newton step in the full $n$-dimensional
space, which is generally too costly when $n$ is huge.

In the general case where $H$ is not necessarily positive definite, an update
of the form

$w\leftarrow w-\eta\,V(V^{T}HV+\alpha I)^{-1}V^{T}g,$ (4)

may be viewed as a regularized Newton step for the restriction of $f$, i.e.,
we find $z^{\star}$ by minimizing $\tilde{h}(z)+\alpha/2\|z\|_{2}^{2}$. This
update is guaranteed to be a descent direction if $V^{T}g\neq 0$ and $\alpha$
is chosen such that $V^{T}HV+\alpha I\succ 0$. We note that

$\displaystyle(\alpha I+VV^{T}HVV^{T})^{-1}$ $\displaystyle=$ $\displaystyle
V(V^{T}HV+$ $\displaystyle\alpha I)^{-1}V^{T}+\alpha^{-1}(I-VV^{T}),$

so if $g\in\operatorname{Range}(V)$, which implies that $(I-VV^{T})g=0$, then
(4) is equivalent to

$w\leftarrow w-\eta\,(\alpha I+VV^{T}HVV^{T})^{-1}g,$

where $VV^{T}HVV^{T}$ is a low-rank approximation to $H$. Moreover, this
update is similar to the update rule used in SketchySGD [6], which is of the
form

$w\leftarrow w-\eta\,(\rho I+\hat{H})^{-1}g,$

where $\hat{H}$ is a low-rank approximation to $H$ that is obtained by means
of a randomized Nyström approximation.

###  II-A Adaptive Subspace Selection

We now consider an iteration of the form

$\displaystyle
w_{k+1}=w_{k}-\eta_{k}\,V_{k}(V_{k}^{T}H_{k}V_{k}+\alpha_{k}I)^{-1}V_{k}^{T}g_{k}$
(5)

for $k=0,1,2,\ldots$ and where $w_{0}\in\mathbb{R}^{n}$ is an initial guess,
$g_{k}=\nabla f(w_{k})$, $H_{k}=\nabla^{2}f(w_{k})$, and
$V_{k}\in\mathbb{R}^{n\times m}$ has orthonormal columns. We will construct
$V_{k}$ such that $g_{k}\in\operatorname{Range}(V_{k})$, and we choose
$\alpha_{k}$ such that $V_{k}^{T}H_{k}V_{k}+\alpha_{k}I\succ 0$. This yields a
descent direction, whenever $\nabla f(w_{k})\neq 0$, since

$\displaystyle\left.\frac{d}{d\eta}f(w_{k}-\eta
V_{k}(V_{k}^{T}H_{k}V_{k}+\alpha_{k}I)^{-1}V_{k}^{T}g_{k})\right|_{\eta=0}$
$\displaystyle=-\|V_{k}^{T}g_{k}\|_{(V_{k}^{T}H_{k}V_{k}+\alpha_{k}I)^{-1}}^{2}<0,$

where $\|x\|^{2}_{W}=x^{T}Wx$ denotes the weighted norm of $x$ induced by the
symmetric positive definite matrix $W$ of order $n$. For example, $V_{k}$ can
be constructed such that

$\operatorname{Range}(V_{k})\supseteq\operatorname{Span}(\nabla
f(w_{k}),\nabla f(w_{k-1}),\dots,\nabla f(w_{k-m+1})),$

and $\alpha_{k}$ can be chosen as
$\alpha_{k}=\max(0,\rho-\lambda_{\min}(V_{k}^{T}H_{k}V))$ where $\rho>0$ is a
given threshold, which implies that
$\lambda_{\min}(V_{k}^{T}H_{k}V_{k}+\alpha_{k}I)\geq\rho$.

###  II-B Block-Diagonal Approximation for Deep Neural Networks

For deep NNs, the computation of Hessian-vector products can be quite
expensive, and to reduce this cost, we may consider an approximation in which
the Hessian is computed layer for layer. This corresponds to a block-diagonal
approximation to the Hessian.

Given a NN with $p$ layers, we partition the parameter vector $w_{k}$ as
$w_{k}=(w_{k,1},\ldots,w_{k,p})$ such that $w_{k,l}$ contains all the
parameters associated with layer $l$. Similarly, if we let
$f\colon\mathbb{R}^{n}\to\mathbb{R}$ denote the loss function and partition
$\nabla f(w_{k})$ as $g_{k}=(g_{k,1},\ldots,g_{k,p})$ where $g_{k,l}$ is the
gradient of $f(w_{k})$ with respect to $w_{k,l}$, then

$H_{k,l}v_{l}=\frac{\partial(g_{k,l}^{T}v_{l})}{\partial w_{k,l}}$

is the product of the Hessian associated with layer $l$ and a vector $v_{l}$.
The resulting approximation to $H_{k}$ can be expressed as

$\hat{H}_{k}=\begin{bmatrix}H_{k,1}&0&\dots&0\\\ 0&H_{k,2}&\dots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\dots&H_{k,p}\\\ \end{bmatrix},$

and a matrix-vector product $\hat{H}_{k}v$ can be computed efficiently using
automatic differentiation. Algorithm 1 outlines a variant of the AdaSub
method, which we will analyze in a stochastic setting in the next section.

Algorithm 1 AdaSub with block-diagonal Hessian approx.

subspace dimension $m$, learning rate $\eta>0$, threshold $\rho>0$

$k\leftarrow 0$

repeat

for $l=1$ to $p$ do

if $k=0$ then

$G_{l}\leftarrow g_{k,l}$

else

$G_{l}\leftarrow\begin{bmatrix}G_{l}&g_{k,l}\end{bmatrix}$

if $k\geq m$ then

Delete the first column of $G_{l}$

end if

end if

Compute orthonormal basis $V_{k,l}$ for $\operatorname{Range}(G_{l})$

Compute $Y\leftarrow H_{k,l}V_{k,l}$

Form $Y^{T}V_{k,l}$

Compute spectral decomposition $Y^{T}V_{k,l}=U\Lambda U^{T}$

$\alpha\leftarrow\max(0,\rho-\min(\Lambda))$

$w_{k+1,l}\leftarrow w_{k,l}-\eta_{k}\,V_{k,l}U(\Lambda+\alpha
I)^{-1}U^{T}V_{k,l}^{T}g_{k,l}$

end for

$k\leftarrow k+1$

until stopping criteria met

##  III Analysis

We now analyze the iteration (5) in a stochastic setting under the assumption
that the function $f$ is twice continuously differentiable and of the form

$\displaystyle f(w)=\frac{1}{N}\sum_{i=1}^{N}f_{i}(w)+\phi(w),$ (6)

where $\phi\colon\mathbb{R}^{n}\to\mathbb{R}$ is a regularization function.
Given two random subsets $S_{k}$ and $\tilde{S}_{k}$ of $\\{1,\ldots,N\\}$, we
define a minibatch gradient and Hessian at iteration $k$ as

$\displaystyle g_{k}$
$\displaystyle=\frac{1}{|S_{k}|}\sum_{i=1}^{|S_{k}|}\nabla
f_{i}(w_{k})+\nabla\phi(w_{k}),$ $\displaystyle H_{k}$
$\displaystyle=\frac{1}{|\tilde{S}_{k}|}\sum_{i=1}^{|\tilde{S}_{k}|}\nabla^{2}f_{i}(w_{k})+\nabla^{2}\phi(w_{k}).$

Furthermore, we will make the following assumptions:

######  Assumption 1.

There exist positive constants $\mu$ and $L$ such that for all
$w\in\mathbb{R}^{n}$,

$\mu I\preceq\nabla^{2}f(w)\preceq LI,$

i.e., $f$ is $\mu$-strongly convex and $L$-smooth.

######  Assumption 2.

The stochastic gradients are unbiased, i.e.,
$\operatorname{\mathbb{E}}[\,g_{k}\,|\,w_{k}]=\nabla f(w_{k})$, the stochastic
Hessians satisfy $\mu I\preceq H_{k}\preceq LI$, and $g_{k}$ and $H_{k}$ are
independent random variables.

######  Assumption 3.

The stochastic gradients are bounded in the sense that [8]

$\operatorname{\mathbb{E}}[\,\|g_{k}\|_{2}^{2}\,|\,w_{k}\,]\leq 2\xi
L(f(w_{k})-f(w^{\star}))+\sigma^{2}$ (7)

where $\sigma^{2}=\operatorname{\mathbb{E}}[\,\|g(w^{\star})\|_{2}^{2}\,]$ is
the gradient noise of $f$ at the optimum and $\xi$ is a positive constant such
that $\operatorname{\mathbb{E}}[\,\|g_{k}\|_{2}^{2}\,]\leq\xi\|\nabla
f(w_{k})\|_{2}^{2}$. The inequality (7) is known as the “weak growth
condition” [15] in the special case where $\sigma^{2}=0$.

Assumption 1 implies that for all $w$ and $w^{\prime}$,

$f(w^{\prime})\leq f(w)+\nabla
f(w)^{T}(w^{\prime}-w)+\frac{L}{2}\|w^{\prime}-w\|_{2}^{2},$ (8)

and strong convexity implies that

$\|\nabla f(w)\|_{2}^{2}\geq 2\mu(f(w)-f(w^{\star})),$ (9)

which is known as the Polyak–Lojaciewicz (PL) inequality [10]. We note that
the independence of $g_{k}$ and $H_{k}$ in Assumption 2 can be ensured by
letting the subsets $S_{k}$ and $\tilde{S}_{k}$ be independent.

######  Theorem 1.

Let $f\colon\mathbb{R}^{n}\to\mathbb{R}$ be a twice continuously
differentiable function, instate Assumptions 1–3, and let $\kappa=L/\mu$. The
iteration (5) with the step size $\eta_{k}=\eta\in(0,\kappa^{-3}\xi^{-1}]$
then satisfies

$\displaystyle\operatorname{\mathbb{E}}[\,f(w_{k})\,]-f(w^{\star})\leq\left(1\\!-\\!\frac{\eta}{\kappa}\right)^{k}(f(w_{0})-f(w^{\star}))+\frac{\kappa^{2}\eta}{2\mu}\sigma^{2}.$
(10)

###### Proof.

From the iteration (3), we have that
$w_{k+1}-w_{k}=-\eta\tilde{H}_{k}^{+}g_{k}$ where
$\tilde{H}_{k}^{+}=V_{k}(V_{k}^{T}H_{k}V_{k})^{-1}V_{k}^{T}$ denotes the
Moore–Penrose pseudo inverse of
$\tilde{H}_{k}=V_{k}(V_{k}^{T}H_{k}V_{k})V_{k}^{T}$. Moreover, using the
inequality (8), we have

$\displaystyle f(w_{k+1})$ $\displaystyle\leq f(w_{k})-\eta\,\nabla
f(w_{k})^{T}\tilde{H}_{k}^{+}g_{k}+\frac{L\eta^{2}}{2}\|\tilde{H}_{k}^{+}g_{k}\|_{2}^{2}$
$\displaystyle\leq f(w_{k})-\eta\,\nabla
f(w_{k})^{T}\tilde{H}_{k}^{+}g_{k}+\frac{L\eta^{2}}{2}\|\tilde{H}_{k}^{+}\|_{2}^{2}\|g_{k}\|_{2}^{2}.$

Assumption 2 implies that
$\|\tilde{H}_{k}^{+}\|_{2}=\|(V_{k}^{T}H_{k}V_{k})^{-1}\|_{2}\leq 1/\mu$,
which leads to

$\displaystyle f(w_{k+1})$ $\displaystyle\leq f(w_{k})-\eta\,\nabla
f(w_{k})^{T}\tilde{H}_{k}^{+}g_{k}+\frac{L\eta^{2}}{2\mu^{2}}\|g_{k}\|_{2}^{2}.$
(11)

Taking the expectation and using Assumptions 2–3, we find that

$\displaystyle\operatorname{\mathbb{E}}[\,f(w_{k+1})\,]$ $\displaystyle\leq
f(w_{k})-\eta\,\nabla
f(w_{k})^{T}\operatorname{\mathbb{E}}[\,\tilde{H}_{k}^{+}\,]\operatorname{\mathbb{E}}[\,g_{k}\,]$
$\displaystyle\qquad+\frac{L\eta^{2}}{2\mu^{2}}\operatorname{\mathbb{E}}[\,\|g_{k}\|_{2}^{2}\,]$
$\displaystyle\leq f(w_{k})-\frac{\eta}{L}\,\|\nabla
f(w_{k})\|_{2}^{2}+\frac{L\eta^{2}}{2\mu^{2}}\operatorname{\mathbb{E}}[\,\|g_{k}\|_{2}^{2}\,],$

where the last inequality follows from Assumption 2 by noting that
$\|\operatorname{\mathbb{E}}[\,\tilde{H}_{k}^{+}\,]\|_{2}\geq 1/L$. Next,
using the PL inequality (9) and Assumption 3, we get

$\displaystyle\operatorname{\mathbb{E}}[\,f(w_{k+1})\,]-f(w^{\star})\leq$
$\displaystyle\qquad\left(1-\frac{2\eta\mu}{L}+\frac{L^{2}\eta^{2}\xi}{\mu^{2}}\right)(f(w_{k})-f(w^{\star}))+\frac{L\eta^{2}}{2\mu^{2}}\sigma^{2}.$

Letting $\eta\leq\mu^{3}/(L^{3}\xi)$, we see that

$1-\frac{2\eta\mu}{L}+\frac{L^{2}\eta^{2}\xi}{\mu^{2}}=1-\frac{2\eta\mu}{L}\left(1-\frac{L^{3}\eta\xi}{2\mu^{3}}\right)\leq
1-\frac{\eta\mu}{L}=1-\frac{\eta}{\kappa},$

and hence,

$\displaystyle\operatorname{\mathbb{E}}[\,f(w_{k+1})\,]-f(w^{\star})\leq\left(1\\!-\\!\frac{\eta}{\kappa}\right)(f(w_{k})-f(w^{\star}))+\frac{\kappa\eta^{2}}{2\mu}\sigma^{2}.$

Finally, we arrive at (10) by recursively taking the expectation over all
steps. ∎

A direct consequence of Theorem 1 is that
$\operatorname{\mathbb{E}}[\,f(w_{k})\,]-f(w^{\star})\leq\epsilon$ provided
that the learning rate $\eta$ is sufficiently small and $k$ is sufficiently
large.

Although analyzing the non-convex case is important, it is a complex problem
that requires further investigation. Thus, we will defer the analysis of the
non-convex case to future work.

##  IV Numerical Experiments

To assess the practical performance of the AdaSub method, we have conducted
four experiments based on different settings, including convex and non-convex
problems and different architectures of NNs. In the first experiment, we apply
the method to a logistic regression problem using a Python implementation of
the method based on the JAX [2] library. The remaining experiments utilize
NNs, for which we implement a PyTorch [14] optimizer that incorporates the
block-diagonal Hessian approximation described in §II-B. We also use the ReLU
activation function, despite its lack of twice continuous differentiability.
In all our experiments, we use $\rho=0.02$ as the lower bound for the
eigenvalues of $V_{k}^{T}H_{k}V_{k}+\alpha_{k}I$.

We compare our method to widely used methods such as Adam [11], AdaGrad [5],
and SGD in all experiments. For the experiments involving a NN, we also
compare to AdaHessian [16]. We note that we were unable to obtain meaningful
results with SketchySDG’s default parameter choices, and therefore, we defer a
fair comparison with SketchySDG to future work.

The learning rates for all the experiments were taken from papers with similar
numerical experiments [16, 6] or manually chosen to ensure the best
performance of each optimizer and a fair comparison.

###  IV-A Logistic Regression

In our first experiment, we trained a logistic regression model to classify if
a given image features a cat. We used a data set with 1200 cat images from
Kaggle [7] and 1200 non-cat images of random objects (such as food, cars,
houses, etc) from the ”other” category of the Linnaeus 5 dataset [3]. Each
image has $64\times 64$ pixels, and we used 80% of the images for training and
20% for testing. We used a batch size of 192 images, $\eta=0.2$, and a
subspace of dimension $m=10$. With Adam/AdaGrad/SGD, the learning rates were
tuned manually and set to be 0.005/0.001/0.2, and with Adam, we used the
parameters $\beta_{1}=0.9$ and $\beta_{2}=0.999$.

The ensure a lower bound for the stochastic Hessian we used a weight decay
with regularization parameter $\lambda=10^{-4}$ in this experiment. Figure 1
shows the results from this experiment and demonstrates AdaSub’s advantage
over the other methods in terms of both training loss and test accuracy.

Figure 1: Accuracy of the logistic regression model in the test set per epoch
and value of the cost function in the training set per training step.

###  IV-B Fully Connected Neural Network

Next, to test our method in a fully connected NN, we built a shallow fully
connected NN and trained it on the MNIST [13] dataset. This dataset consists
of 70,000 grayscale images of handwritten numbers. We used 60,000 images as
the training set and the remaining 10,000 images as the test set, and we used
a batch size of 1000 images. Each image is represented as a vector of length
784 (corresponding to 28x28 pixels). We used a network with a single hidden
layer with 100 units and a ReLU activation function, and as the output layer,
we used a softmax activation function with 10 outputs. For AdaSub, we used
$\eta=0.1$ and a subspace of dimension $m=40$. For
AdaHessian/Adam/AdaGrad/SGD, the learning rates were tuned manually and set to
be 0.05/0.001/0.01/0.1. For AdaHessian/Adam, we used $\beta_{1}=0.9$ and
$\beta_{2}=0.999$. We note that this is a nonconvex problem, and hence the
analysis from the previous section does not apply.

Figure 2 shows, for all the methods, the accuracy in the test set after each
epoch and the cost function for the current minibatch at each iteration.
AdaSub presented an advantage over the other methods in terms of both training
loss and test accuracy.

Figure 2: Accuracy of the fully connected network in the test set per epoch
and value of the cost function in the training set per training step.

###  IV-C Convolutional Neural Network

We have also applied our method to the problem of training a convolutional
neural network (CNN) based on the CIFAR-10 [12] dataset. This dataset consists
of 60,000 32x32 color images divided into 10 different classes. The classes
are airplanes, cars, birds, cats, deer, dogs, frogs, horses, ships, and
trucks, and each class has the same number of examples. We used 50,000 images
as the training set, 10,000 images as the test set, and we used a batch size
of 400 images.

The input to the CNN is a color image with dimensions 32x32x3. The network
begins with a convolutional layer that has six output channels and uses 5x5x3
dimensional kernels. After the convolutional layer, a 2x2 max pooling layer
with a stride of 2 is applied. This process is repeated with a second
convolutional layer that has 16 output channels and uses 5x5x6 dimensional
kernels, followed by another 2x2 max pooling layer with a stride of 2. The
output from this step is a tensor with dimensions of 5x5x16. This tensor is
flattened into a vector of length 400, which is then fed into a fully
connected network with two hidden layers. The first hidden layer has 120
units, and the second has 84 units. The output layer has 10 units and uses a
softmax activation function. All the layers except the output layer use a ReLU
activation function.

In this experiment, we used $\eta=0.15$ and a subspace of dimension $m=20$ in
our method. For AdaHessian/Adam/AdaGrad/SGD, the learning rates were tuned and
set to be 0.1/0.001/0.01/0.4. For AdaHessian/Adam, we used $\beta_{1}=0.9$ and
$\beta_{2}=0.999$.

Figure 3 shows, for all the methods, the accuracy in the test set after each
epoch and the cost function for the current minibatch at each iteration.
Comparing these results with the ones from the previous experiment, we can
observe that even for a more complex problem with a smaller adaptive subspace,
AdaSub was able to outperform the other methods in terms of both training loss
and test accuracy.

Figure 3: Accuracy of the CNN in the test set per epoch and value of the cost
function in the training set per training step.

###  IV-D ResNet20

To test our method for the training of deeper NNs, we applied it to ResNet20
[9] using the CIFAR-10 dataset. We used $\eta=0.15$ and a batch size of 400
images. Since ResNet20 is a NN with approximately 270,000 parameters, we chose
to use a two-dimensional subspace ($m=2$) for this experiment to limit the
computational cost. Thus, we only need to perform two (approximate) Hessian-
vector products at each iteration. For AdaHessian/Adam/AdaGrad/SGD, the
learning rates were tuned and set to 0.2/0.001/0.1/0.1. For AdaHessian/Adam,
we used $\beta_{1}=0.9$ and $\beta_{2}=0.999$.

Figure 4 shows, for all the methods, the accuracy in the test set and the cost
function for the current minibatch after each epoch. Comparing these results
with the last experiments, we observe that our method still outperforms
SDG/Adam/Adagrad and is comparable to AdaHessian despite the small subspace
dimension.

Figure 4: Accuracy of ResNet20 in the test set and value of the cost function
in the training set per epoch.

###  IV-E Time Comparison

The cost of one minibatch iteration using AdaSub or AdaHessian is higher than
that of Adam, Adagrad, and SGD because of the added cost of Hessian-vector
products. Thus, to compare the runtime, we repeated the first experiment
(logistic regression), letting each optimizer run for 10 seconds.

Figure 5 shows the test accuracy and the training loss as a function of the
runtime. From the test accuracy plot, we can see that the test accuracy
obtained after 10 seconds using Adam is matched by AdaSub after approximately
1.4 seconds. Moreover, in terms of training loss, AdaSub outperforms Adam
after approximately 1 second.

Figure 5: Test accuracy and training loss as a function of time.

##  V Conclusions and Future Work

We have presented AdaSub, a new stochastic optimization method that restricts
the search direction at each iteration to a low-dimensional subspace of the
previous gradients. Our method achieves a better balance between accuracy and
computational efficiency by incorporating second-order information into the
optimization process while limiting the dimension of the subspace. We tested
the method on both convex and non-convex problems, and in both settings, it
outperformed several state-of-the-art optimization methods. These results
suggest that our proposed algorithm can be a promising approach for optimizing
complex and high-dimensional models in a range of applications within machine
learning. Moreover, our preliminary experiments suggest that AdaSub is also
relatively robust to the choice of algorithmic parameters such as the number
of past directions and the step size.

The current version of the algorithm is a first attempt at incorporating
second-order information into the optimization process without significantly
increasing the cost per iteration. The construction of the subspace and the
practical implications of its dimension need to be further explored. For
example, the construction of $V_{k}$ may be based on additional information
such as previous search directions in addition to previous gradients.
Moreover, a more thorough analysis of the convergence rate of the algorithm is
needed under less restrictive and more realistic assumptions, and additional
numerical experiments are needed to test the performance of the algorithm on a
wider range of problems. In particular, we would like to test the algorithm on
larger datasets and deeper NNs.

## References

  * [1] Albert S. Berahas, Majid Jahani, Peter Richtárik, and Martin Takáč.  Quasi-Newton methods for machine learning: Forget the past, just sample, 2019.
  * [2] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2018.
  * [3] Giorgi Chaladze and Levan Kalatozishvili.  Linnaeus 5 dataset for machine learning, 2017.
  * [4] Damek Davis, Dmitriy Drusvyatskiy, Sham Kakade, and Jason D. Lee.  Stochastic subgradient method converges on tame functions, 2018.
  * [5] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of Machine Learning Research, 12(61):2121–2159, 2011.
  * [6] Zachary Frangella, Pratik Rathore, Shipu Zhao, and Madeleine Udell.  SketchySGD: Reliable stochastic optimization via robust curvature estimates, 2022.
  * [7] Spandan Ghosh.  Cats faces 64x64 (for generative models), Feb 2019.
  * [8] Robert Mansel Gower, Nicolas Loizou, Xun Qian, Alibek Sailanbayev, Egor Shulgin, and Peter Richtárik.  SGD: General analysis and improved rates.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 5200–5209. PMLR, 09–15 Jun 2019.
  * [9] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [10] Hamed Karimi, Julie Nutini, and Mark Schmidt.  Linear convergence of gradient and proximal-gradient methods under the polyak-Łojasiewicz condition, 2016.
  * [11] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2014.
  * [12] Alex Krizhevsky and Geoffrey Hinton.  Learning multiple layers of features from tiny images.  Technical Report 0, University of Toronto, Toronto, Ontario, 2009.
  * [13] Yann LeCun and Corinna Cortes.  MNIST handwritten digit database.  2010\.
  * [14] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [15] Sharan Vaswani, Francis R. Bach, and Mark Schmidt.  Fast and faster convergence of SGD for over-parameterized models and an accelerated perceptron.  CoRR, abs/1810.07288, 2018.
  * [16] Zhewei Yao, Amir Gholami, Sheng Shen, Kurt Keutzer, and Michael W. Mahoney.  ADAHESSIAN: an adaptive second order optimizer for machine learning.  CoRR, abs/2006.00719, 2020.


# WELL-POSEDNESS FOR CAUCHY FRACTIONAL PROBLEM INVOLVING DISCRETE CONVOLUTION
OPERATORS

Jorge González-Camus  Departamento de Matemática, Facultad de Ciencias
Naturales, Matemática y del Medio Ambiente, Universidad Tecnológica
Metropolitana, Santiago, Chile<EMAIL_ADDRESS>
###### Abstract.

This work is focused in establishing sufficient conditions to guarantee the
well-posedness of the following nonlinear fractional semidiscrete model

$\begin{cases}\mathbb{D}^{\beta}_{t}u(n,t)=Bu(n,t)+f(n-ct,u(n,t)),\,&n\in\mathbb{Z},\;t>0,\\\
u(n,0)=\varphi(n),\;&n\in\mathbb{Z},\\\ \end{cases}$

under the assumptions that $\beta\in(0,1]$, $c>0$ some constant, $B$ is a
discrete convolution operator with kernel $b\in\ell^{1}(\mathbb{Z})$, which is
the infinitesimal generator of the Markovian $C_{0}$-semigroup and suitable
nonlinearity $f$. We present results concerning the existence and uniqueness
of solution, as well as establishing a comparison principle of solutions
according to respective initial values and applications.

###### Key words and phrases:

Caputo fractional derivative; discrete fractional Laplacian; discrete
fractional operators; fundamental solutions; Wright and Mittag-Leffler
functions; discrete convolution operator, Markovian Semigroup

###### 2010 Mathematics Subject Classification:

35R11, 35A08, 39A12, 47D07

##  1\. Introduction

Fractional evolution equations have become an important focus of study in the
last decades, due to the numerous phenomena modeled in which the memory effect
is observable. We can found examples of this fact arising from diverse areas
of knowledge, such as physics, engineering, biology, among others. At the same
time, semidiscrete models have received special attention since, from
numerical analysis, it is well known that semidiscretizations is a useful
method to study the general behavior of evolutionary partial differential
equations. This topic has been studied by the branch of mathematical analysis,
using techniques from PDEs, functional and harmonic analysis, among others. We
refer [8, 11, 12, 18, 19, 20, 22, 36, 37, 38] and references therein.

In regards to fractional powers of difference operators, the _discrete
fractional Laplacian_ $(-\Delta_{d})^{\alpha}$ and _fractional backward
/forward Euler operator_ $\Delta^{\alpha}/(-\nabla)^{\alpha}$ have received
special attention, due to the numerous applications, mainly diffusive and
transport nature. We refer recent references [2, 9, 11, 12, 20, 23].

First, we consider the following general linear semidiscrete equation

(1.1)
$\begin{cases}\mathbb{D}^{\beta}_{t}u(n,t)=Bu(n,t)+f(n,t),&n\in\mathbb{Z},\;t>0,\\\
u(n,0)=\varphi(n),\;&n\in\mathbb{Z},\\\ \end{cases}$

where $Bf=b\ast f$ is a discrete convolution operator over
$\ell^{p}(\mathbb{Z}),\ p\in[1,\infty]$. We remark that in the particular case
$B=(-\Delta_{d})^{\alpha}$, we recover several semidiscrete models:

  * •

For $\alpha=\beta=1$ we get the _semidiscrete heat equation_ (see [9]).

  * •

For $\alpha=\frac{1}{2},\beta=1$ we get the _semidiscrete Poisson equation_
(see [9]).

  * •

For $\alpha=1,\beta=\frac{1}{3}$ we get the _semidiscrete Airy equation_ (see
[39]).

Similarly, for $B=(-\nabla)^{\alpha}/\Delta^{\alpha}$, we get the linear
_fractional semidiscrete transport equation_ (see [2]):

(1.2)
$\begin{cases}\mathbb{D}^{\beta}_{t}u(n,t)=-((-\nabla)^{\alpha}/\Delta^{\alpha})u(n,t)+f(n,t),\,&n\in\mathbb{Z},\;t>0,\\\
u(n,0)=\varphi(n),\;&n\in\mathbb{Z}.\\\ \end{cases}$

The study of semidiscrete models is not new. Under our knowledge, the
pioneering work on this topic was presented by Bateman in [5] in 1943, where
the author considered a broad set of discrete operators and the influence on
its $C_{0}$-semigroup associated, including the entire cases for the systems
(1.1) and (1.2), with application to a wide array of physical and engineering
systems. The list presented shows us the crucial role played by the special
functions, such as Bessel type, Hermite’s polynomial, among others. Later,
many researchers focused their efforts on the study of properties of the
difference operators. For instance, for diffusive and transport models, we
refer [1, 32, 33, 34, 35, 40].

Interest about fundamental solutions from models involving fractional powers
of the Laplacian comes partially from an equation studied by Mainardi, Luchko
and Pagnini [24]. Since then, the study of explicit representation of solution
for Cauchy problem involving this and other operators has been focus of
attention of several researchers in the last decade.

In regards to fractional discrete models, Keyantuo et. al. in [20] determined
the explicit representation of fundamental solution for (1.1) in the
subdiffusive fractional case (i.e, $0<\alpha,\beta\leq 1$) through the
subordination principles on Wright and Lévy functions (see [3]), obtaining an
integral representation. Furthermore, the authors proved the fact that the
discrete fractional Laplacian generates a Markovian semigroup (see [20,
Theorem 3.8]). After, González-Camus et. al. in [11] got the representation of
the explicit fundamental solution for (1.1) in the subdiffusive and
superdiffusive fractional cases (i.e, $0<\alpha\leq 1$, with $0<\beta\leq 1$
and $1<\beta\leq 2$ respectively) in terms of the generalized Wright
functions. In both cases, the discrete Fourier transform was the fundamental
tool in order to get their respective representations.

For the fractional backward/forward Euler operator, in [2] the authors
obtained maximum and comparison principles results for the discrete fractional
derivatives, as well as regularity results when the space is a mesh of length
$h$, and approximation results to the continuous fractional derivative in
Marchaud sense.

A general method for the study of the system (1.1) was presented by González-
Camus et. al in [12], where the solution (unique) was obtained explicitly in
terms of Mittag-Leffler on the Banach Algebra $\ell^{1}(\mathbb{Z})$, as well
as an integral representation, based in subordination principles and
functional calculus presented in [3] and [6], respectively. The explicit
representation through Mittag-Leffler function is given by

$\displaystyle\displaystyle u(n,t)=$
$\displaystyle(E_{\beta,1}(t^{\beta}b)\ast\varphi)(n)$
$\displaystyle+\int_{0}^{t}(t-s)^{\beta-1}\left(E_{\beta,\beta}((t-s)^{\beta}b)\ast
f(\cdot,s)\right)(n)\ ds,\,\,n\in\mathbb{Z}.$

We remark that this result generalizes [19, Theorem 3.3] and the case
$0<\alpha\leq 1$ studied by González et. al. in [11, Theorem 5.1].
Furthermore, the authors proved that if the initial data
$\varphi\in\ell^{p}(\mathbb{Z})$, for $1<p\leq\infty$, then this solution $u$
belongs to $\ell^{p}(\mathbb{Z})$.

Nonlinear models, and particularly nonlinear discrete systems, have been the
subject of research of many authors, due to their relevance in such areas as
spatial ecology, physics, chemistry, biology and the emerging area of global
warming. For instance, the Fisher equation

$u_{t}(n,t)=\Delta_{d}u(n,t)+ru(n,t)-u^{2}(n,t),\ n\in\mathbb{Z},\ t\geq 0,$

was first studied by Zinner, Harris and Hudson [44] in the context of
traveling waves. The version with a retarded term, i.e.

$u_{t}(n,t)=\Delta_{d}u(n,t)+r(n-ct)u(n,t)-u^{2}(n,t),\ n\in\mathbb{Z},\ t\geq
0,$

for $c>0$ was studied by Hu and Li in [16]. Another interesting example is
provided by the discrete _Nagumo equation_ , given by

$u_{t}(n,t)=\rho\Delta_{d}u(n,t)+u(n,t)(u(n,t)-a)(1-u(n,t)),\ n\in\mathbb{Z},\
t\geq 0,$

where $\rho>0$ and $a\in(0,\frac{1}{2})$. The discrete Nagumo equation has
been focus of attention due to its differences with the continuous
counterpart. For more information we refer [17, 20].

Due to the above examples of nonlinear models, Keyantuo et. al. in [20]
studied the model

$\begin{cases}\mathbb{D}^{\beta}_{t}u(n,t)=-(\Delta_{d})^{\alpha}u(n,t)+f(n-ct,u(n,t)),&n\in\mathbb{Z},\;t>0,\\\
u(n,0)=\varphi(n),\;&n\in\mathbb{Z},\\\ \end{cases}$

considering a suitable function $f$ and $0<\alpha\leq 1$, obtaining general
results including all equations presented above.

We can note that all the above nonlinear models are governed by the spacial
operator discrete fractional Laplacian. Here arises a natural question: it is
an intrinsic property of the discrete fractional Laplacian, or is there a
deeper concept which is crucial in order to obtain this well-posedness? The
main purpose of this study is to answer the question above. Our purpose in
this paper is to provide sufficient conditions for the well-posedness of the
nonlinear model

$\begin{cases}\mathbb{D}^{\beta}_{t}u(n,t)=Bu(n,t)+f(n-ct,u(n,t)),&n\in\mathbb{Z},\;t>0,\\\
u(n,0)=\varphi(n),\;&n\in\mathbb{Z},\\\ \end{cases}$

where we consider $\beta\in(0,1]$ is real number, $c>0$ some positive
constant, the operator $Bh(n)=(b\ast h)(n)$, with $b$ belonging to Banach
Algebra $\ell^{1}(\mathbb{Z})$ being a infinitesimal generator of a Markovian
semigroup and $f$ and $\varphi$ suitable functions, $p\in[1,\infty]$. We
recall that $\mathbb{D}_{t}^{\beta}$ denotes the _Caputo fractional
derivative_ given by

$\mathbb{D}_{t}^{\beta}v(t)=\frac{1}{\Gamma(1-\beta)}\int_{0}^{t}(t-s)^{-\beta}v^{\prime}(s)ds,$

for $t>0$ and $\alpha>0$, and $v$ a differentiable function, with the object
of extending the result obtained in [20] for a branch of important operators
in differences and its fractional powers. To answer this question, we will use
tools from Banach Algebra $\ell^{1}(\mathbb{Z})$, resolvent families of
operators theory and properties of Markov operators.

This paper is organized as follows: In Section 2, we present the basic
concepts of this work, focused mainly in properties of discrete convolution
operators, the Banach Algebra $\ell^{1}({\mathbb{Z}})$, Markovian semigroups
and resolvent operators theory. In Section 3 we present the main result of
this study in the Theorem (3.3). We give the sufficient conditions to ensure
the well-posedness of the nonlinear model (3.1) and uniqueness of solution,
through an iterative method and Gronwall inequality. Finally, we present two
corollaries concerning the comparison principle for solutions according to the
respective initial values. In Section 4 we give some applications for the main
result obtained in the Section 3 and get the Markovian property for semigroup
generated by fractional powers of a convolution operator which is an
infinitesimal generator of a Markovian semigroup, presented in the Theorem
(4.1). We finish this paper with an Appendix, with the definitions and
principal properties about special functions used throughout this study.

##  2\. Preliminaries

In this section, we give the basic concepts and results concerning to discrete
convolution operator. Furthermore, we present the notation used throughout
this work.

###  2.1. Operators defined by discrete convolution

We consider as framework of this study the commutative _Banach Algebra with
unit_

$\ell^{1}(\mathbb{Z}):=\\{f:\mathbb{Z}\to\mathbb{R}:\sum_{n\in\mathbb{Z}}|f(n)|<\infty\\},$

equipped with the product the _discrete convolution_ defined by

$(f\ast g)(n):=\sum_{j\in\mathbb{Z}}f(n-j)g(j),\qquad n\in\mathbb{Z}.$

The unit for the product is given by _$\delta$ -Kronecker measure_
concentrated at $0$, denoted by $\delta_{0}$. This element satifies: for all
$a\in\ell^{1}(\mathbb{Z}),\ a\ast\delta_{0}=a$. It is convenient recall that
the $\delta$-Kronecker measure $\delta_{n}$ is defined by

$\delta_{n}(m)\begin{cases}1,&\quad m=n,\\\ 0,&\quad m\neq n.\end{cases}$

Furthermore, we remark that $\ell^{1}(\mathbb{Z})$ is semi-simple and regular
Algebra (see [21, Theorem 4.7.4]) and [21, Corolary 7.2.3] respectively).

Given $a=(a(n))_{n\in\mathbb{Z}}\in\ell^{1}(\mathbb{Z})$, we define the
following discrete operator by convolution

$Af(n):=(a\ast f)(n),\ {n\in\mathbb{Z}},\ f\in\ell^{p}(\mathbb{Z}).$

We say that $A$ is a _discrete convolution operator_ , with _kernel_ $a$. We
note that
$A\in{\mathcal{B}}(\ell^{p}(\mathbb{Z})):=\\{B:\ell^{p}(\mathbb{Z})\to\ell^{p}(\mathbb{Z}):B\
\mbox{is linear and bounded}\\}$ for all $1\leq p\leq\infty$, since from
Young’s inequality

$\|Af(n)\|_{p}=\|(a\ast f)(n)\|_{p}\leq\|a\|_{1}\|f\|_{p}.$

Moreover, the operator $A$ satisfies

(2.1) $\|A\|_{1}=\|a\|_{1}.$

Here, we understand the powers for $a\in\ell^{1}(\mathbb{Z})$ through discrete
convolution, given by

$a^{j}:=a^{j}(n):=(\underbrace{a\ast...\ast a}_{j-times})(n).$

Obviously, we have $a^{0}=\delta_{0}$, $a^{1}=a$.

In the Banach Algebra $\ell^{1}(\mathbb{Z})$, we define $e^{at}$ _the
exponential of $a$_ as follows

$e^{ta}:=\sum_{n=0}^{\infty}{t^{n}a^{n}\over n!},\ t\in\mathbb{R}_{0}^{+}.$

We denote $\\{e^{tA}\\}_{t\geq 0}$ to the _family semigroup of operators_
generated by the operator $A\in\mathcal{B}(\ell^{p}(\mathbb{Z}))$. In the case
of discrete convolution operators, $Af=a\ast f$, this family acts in fact by
discrete convolution over $\ell^{p}(\mathbb{Z})$. Moreover, we have the
explicit representation

$e^{tA}f(n)=(e^{at}\ast f)(n),\ n\in\mathbb{Z},\ f\in\ell^{p}(\mathbb{Z}),$

and, for each $t\geq 0$, from the equality (2.1) we have

(2.2) $\|e^{tA}\|_{1}=\|e^{ta}\|_{1}.$

See [12] and [21] for more information about the Banach Algebra
$\ell^{1}(\mathbb{Z})$ and result for discrete convolution operators. We refer
[10, 27] for more information about the semigroup theory.

Next, we present three important cases that will be used as applications of
this study.

####  2.1.1. The discrete Laplacian

Denoted by $\Delta_{d}$ given by

$\Delta_{d}f(n):=f(n+1)-2f(n)+f(n-1)=([\delta_{-1}-2\delta_{0}+\delta_{1}]\ast
f)(n).$

For $\alpha\in(0,1)$ the _discrete fractional Laplacian_ is given by

$\displaystyle(-\Delta_{d})^{\alpha}f(n)=\sum_{k\in\mathbb{Z}}K_{d}^{\alpha}(n-k)f(k)=(K_{d}^{\alpha}\ast
f)(n),\,n\in\mathbb{Z},$

where the coefficients $K_{d}^{\alpha}$ are given by

$K_{d}^{\alpha}(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}(4\sin^{2}(\theta/2))^{\alpha}e^{-in\theta}d\theta=\frac{(-1)^{n}\Gamma(2\alpha+1)}{\Gamma(1+\alpha+n)\Gamma(1+\alpha-n)},\,n\in\mathbb{Z}.$

We mention that this operator corresponds to fractional powers of the second-
order central difference approximation for the second-order derivative. For
more information about the representation and properties of the discrete
fractional Laplacian, see [11], [12] and [20].

####  2.1.2. The forward/backward Euler operators

Denoted by $-\Delta$ and $\nabla$ respectively, are given by

$-\Delta f(n):=f(n)-f(n+1)=([\delta_{0}-\delta_{-1}]\ast f)(n).$ $\nabla
f(n):=f(n)-f(n-1)=([\delta_{0}-\delta_{1}]\ast f)(n).$

From the above definition, $\alpha\in(0,1)$ the _fractional forward/backward
Euler operator_ , given by

$\displaystyle(\Delta/-\nabla)^{\alpha}f(n)=\sum_{k\in\mathbb{Z}}K_{\pm}^{\alpha}(n-k)f(k)=(K_{\pm}^{\alpha}\ast
f)(n),\,n\in\mathbb{Z},$

where the coefficients $K_{\pm}^{\alpha}$ are given by

$K_{\pm}^{\alpha}(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}(1-e^{\pm
i\theta})^{\alpha}e^{-in\theta}d\theta=(-1)^{n}\frac{\Gamma(\alpha+1)}{\Gamma(\alpha\mp
n+1)\Gamma(\pm n+1)}\chi_{\mathbb{Z}_{\pm}}(n),$

where $n\in\mathbb{Z}$ and $\mathbb{Z}_{\pm}:=\\{\ldots,\pm 2,\pm 1,0\\}$. It
is convenient remark that the kernels $K_{+}^{\alpha}$ and $K_{-}^{\alpha}$
correspond to $\Delta^{\alpha}$ and $(-\nabla)^{\alpha}$ respectively. We
recall that the above operators are related to fractional powers of the Euler
scheme of approximation. For more information about the representation and
properties of the discrete fractional differences see [2] and [12].

####  2.1.3. The discrete Laplacian 2-step

Denoted by $\Delta_{dd}$ given by

$\Delta_{dd}f(n):=f(n+2)-2f(n)+f(n-2)=([\delta_{-2}-2\delta_{0}+\delta_{2}]\ast
f)(n).$

For $\alpha\in(0,1)$ we obtain the _discrete fractional Laplacian 2-step_ ,
given by

$\displaystyle(-\Delta_{dd})^{\alpha}f(n)=\sum_{k\in\mathbb{Z}}K_{dd}^{\alpha}(n-k)f(k)=(K_{dd}^{\alpha}\ast
f)(n),\,n\in\mathbb{Z},$

where the coefficients $K_{dd}^{\alpha}$ are given by

$K_{dd}^{\alpha}(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}(4\sin^{2}(\theta))^{\alpha}e^{-in\theta}d\theta=\frac{\cos({n\over
2}\pi)\Gamma(2\alpha+1)}{\Gamma(1+\alpha+{n\over 2})\Gamma(1+\alpha-{n\over
2})},$

with $n\in\mathbb{Z}$. This operator appears in De Juhasz equation, appearing
in the seminal Bateman’s paper [5] in connection with surges in springs and
connected systems of springs. For more information about the representation
and properties of the discrete fractional Laplacian 2 step, see [12].

The remarkable fact is, for $\alpha\in(0,1)$, $(-\Delta_{d})^{\alpha}$,
$(\Delta/-\nabla)^{\alpha}$ and $(-\Delta_{dd})^{\alpha}$ are linear bounded
operators over $\ell^{p}(\mathbb{Z}),0\leq p\leq\infty$, since for
$|n|\to\infty$, $|K_{d}^{\alpha}(n)|\sim\frac{1}{n^{2\alpha+1}}$,
$|K_{\pm}^{\alpha}(n)|\sim\frac{1}{n^{\alpha+1}}$ and
$|K_{dd}^{\alpha}(n)|\sim\frac{1}{n^{2\alpha+1}}$ respectively. See [12,
Section 4], [25, eq. 4.5, section 4.2] for more information.

It is well-known that the semigroups generated by $\Delta_{d}$,
$-\Delta/\nabla$ and $\Delta_{dd}$ are (see [5], [12])

$e^{t\Delta_{d}}f(n)=\sum_{m\in\mathbb{Z}}e^{-2t}I_{n-m}(2t)f(m)=:(e^{t[\delta_{-1}-2\delta_{0}+\delta_{1}]}\ast
f)(n).$
$e^{t(-\Delta/\nabla)}f(n)=\sum_{m\in\mathbb{Z}}e^{-t}\frac{t^{\mp(n-m)}}{(\mp(n-m))!}\chi_{\mathbb{Z}_{\mp}}(n-m)f(m)=:(e^{t([\delta_{0}-\delta_{-1}]/[\delta_{0}-\delta_{1}])}\ast
f)(n).$
$e^{t\Delta_{dd}}f(n)=\sum_{m\in\mathbb{Z}}e^{-2t}I_{\frac{n-m}{2}}(2t)\chi_{2\mathbb{Z}}(n-m)f(m)=:(e^{t[\delta_{-2}-2\delta_{0}+\delta_{2}]}\ast
f)(n).$

######  Remark 2.1.

The fractional powers for the discrete operators $\Delta_{d}$,
$(-\Delta/\nabla)$ and $\Delta_{dd}$ have been obtained using the
_Balakrishnan’s Formula_. The work [12] a general representation for
fractional powers of discrete operators which are generators of uniformly
bounded $C_{0}$-semigroup in the Banach Algebra $\ell^{1}(\mathbb{Z})$
framework is presented.

###  2.2. Markovian semigroups

In this subsection we present the basic results of a Markovian semigroups used
in this work.

######  Definition 2.2.

Let $(X,\Sigma,\mu)$ be a $\sigma$-finite measure space. We denote
$M\subseteq\ell^{1}(X):=\ell^{1}(X,\Sigma,\mu)$ as the subset of all density
functions over $X$

$M=\\{f\in\ell^{1}(X):f\geq 0,\|f\|_{1}=1\\}.$

We call _Markov Operator_ to all linear mapping $F:\ell^{1}(X)\to\ell^{1}(X)$
satisfying $F(M)\subset M$.

######  Definition 2.3.

A family $\\{T(t)\\}_{t\geq 0}$ of Markov operators is called _Markovian
semigroup_ if furthermore it is a $C_{0}$-semigroup family of operators.

######  Remark 2.4.

We note that if $A$ is a convolution operator with kernel
$a\in\ell^{1}(\mathbb{Z})$, with $e^{tA}$ a Markovian semigroup, then $e^{ta}$
satisfies $e^{ta}\geq 0$ and $\|e^{ta}\|_{1}=1$. Reciprocally, the assertion
is true using (2.2).

######  Remark 2.5.

Let $\textbf{1}(n)$ the sequence constantly equals to $1$. Let $A$ be a
convolution operator with kernel $a\in\ell^{1}(Z)$ be a generator of a
Markovian semigroup $e^{tA}$. Then the following equality holds

(2.3) $e^{tA}\textbf{1}(n)=\textbf{1}(n).$

###### Proof.

We remark that $\textbf{1}(n)=\textbf{1}(n-m)$, for all $n,m\in\mathbb{Z}$.
Since $e^{tA}$ is a Markovian semigroup, we have

$\displaystyle
e^{tA}\textbf{1}(n)=(e^{ta}\ast\textbf{1})(n)=\sum_{m\in\mathbb{Z}}e^{ta}(m)\textbf{1}(n-m)=\textbf{1}(n)\sum_{m\in\mathbb{Z}}e^{ta}(m)=\textbf{1}(n).$

∎

###  2.3. Resolvent and Integral Resolvent families of Operators

A powerful tool in order to study existence of solution of abstract evolution
equations is the _resolvent family of operators theory_ , playing a crucial
role in the study of the well-posedness of fractional systems (see [6, 18,
29]).

######  Definition 2.6.

Let $B$ be a closed linear operator with domain $D(B)\subset X$ and $a\in
L^{1}_{loc}(\mathbb{R}^{+}_{0}).$ A family $\\{S(t)\\}_{t\geq 0}$ of bounded
and linear operators in $X$ is called a _resolvent_ with generator $B$ if the
following conditions are satisfied:

  1. (1)

$S(t)$ is strongly continuous on $\mathbb{R}_{0}^{+}$ and $S(0)=I$;

  2. (2)

$S(t)$ commutes with $B$. That means that $S(t){D}(B)\subset{D}(B)$ and
$BS(t)x=S(t)Bx$, for all $x\in{D}(B)$ and $t\geq 0$;

  3. (3)

The resolvent equation holds:

$S(t)x=x+\int_{0}^{t}a(t-s)BS(s)xds,\quad\mbox{for all}\ x\in{D}(B),\,\,t\geq
0.$

######  Definition 2.7.

Let $B$ be a closed linear operator with domain $D(B)\subset X$ and $a\in
C(\mathbb{R}^{+}_{0}).$ A strongly continuous family $\\{P(t)\\}_{t\geq 0}$ of
bounded linear operators in $X$ is called an _integral resolvent_ with
generator $B$ if the following conditions are satisfied:

  1. (1)

$P(0)=a(0)I;$

  2. (2)

$P(t)$ commutes with $B;$

  3. (3)

The integral resolvent equation holds:

$P(t)x=a(t)x+\int_{0}^{t}a(t-s)BP(s)xds,\quad\mbox{for all}\ x\in{D}(B),t\geq
0.$

######  Remark 2.8.

It is well known the fact for the case $a(t)\equiv 1$ the resolvent family
corresponds to the $C_{0}-$semigroup generated by $A$ (see [29]).

##  3\. Main Results

In this section we present the main purpose of this paper, which is to
establish the well-posedness of the nonlinear model

(3.1)
$\begin{cases}\mathbb{D}^{\beta}_{t}u(n,t)=Bu(n,t)+f(n-ct,u(n,t)),&n\in\mathbb{Z},\;t>0,\\\
u(n,0)=\varphi(n),\;&n\in\mathbb{Z},\\\ \end{cases}$

where $\mathbb{D}_{t}^{\beta}$ is the fractional Caputo deritative,
$\beta\in(0,1]$, $c>0$ some positive constant, the discrete convolution
operator $Bh(n)=(b\ast h)(n)$, with $b$ belonging to Banach Algebra
$\ell^{1}(\mathbb{Z})$ and $f,\varphi$ suitable functions.

In order to prove the main result of this work, for each $\gamma>0$, we define
the set

$\mathcal{L}_{\gamma}:=\\{u\in\ell^{\infty}(\mathbb{Z}):0\leq
u(n)\leq\gamma,\forall n\in\mathbb{Z}\\}.$

Let $\rho>0$. We define the function:
$F:\mathbb{R}\times[0,\gamma]\to\mathbb{R}$, given by

(3.2) $F(x,s):=\rho s+f(x,s).$

This function satisfies properties presented in the following Lemma. The proof
can be found in [20, Lemma 4.1].

######  Lemma 3.1.

Let $\gamma>0$. Suppose that $f:\mathbb{R}\times\mathbb{R}_{0}^{+}$ is
measurable, the restriction of $f(x,\cdot)$ to $[0,\gamma]$ belongs to
$C^{1}[0,\gamma]$, is concave in $[0,\gamma]$, and satisfies $f(x,0)=0$, for
all $x\in\mathbb{R}$. Let $\rho>0$ such that

(3.3) $\rho+\partial_{s}f(x,\gamma)\geq 0,\ \forall x\in\mathbb{R}.$

Then, the function $F$ defined in (3.2) is no decreasing in the second
variable and no negative for $x\in\mathbb{R}$.

We get the expression of the solution for the equation (3.1) employing the
resolvent operator theory and the subordination principle on Wright function
(5.1) presented in [6, 18]. Let us define the following operators

(3.4) $S_{\beta}(t):=\int_{0}^{t}\Phi_{\beta}(\tau)e^{\tau t^{\beta}B}d\tau,\
t\geq 0,$

called the _resolvent family_ for (3.1), and

(3.5) $P_{\beta}(t):=\int_{0}^{t}\beta\tau\Phi_{\beta}(\tau)e^{\tau
t^{\beta}B}d\tau,\ t\geq 0,$

called the _integral resolvent family_ for (3.1). Then, the mild solution
$u(n,t)$ for the nonlinear model (3.1) satisfies the equation

(3.6) $u(n,t)=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)f(n-cs,u(n,s))\
ds.$

We consider the perturbed operator $B_{\rho}:=B-\rho I$. This operator is
bounded over $\ell^{p}(\mathbb{Z})$. Let us denote $\\{e^{tB_{\rho}}\\}_{t\geq
0}$ the semigroup generated by $B_{\rho}$. From general theory of
$C_{0}$-semigroup family of operators, we have $e^{tB_{\rho}}=e^{t(B-\rho
I)}=e^{tB}e^{-\rho t}$. A change of variable of (3.4) and (3.5) gives us an
equivalent representation, given by

(3.7) $S_{\beta}(t):=\int_{0}^{t}t^{-\beta}\Phi_{\beta}(\tau
t^{-\beta})e^{\tau B_{\rho}}d\tau,\ t\geq 0,$

and

(3.8)
$P_{\beta}(t):=\int_{0}^{t}\beta\frac{\tau}{t^{1+\beta}}\Phi_{\beta}(\tau
t^{-\beta})e^{\tau B_{\rho}}d\tau,\ t\geq 0.$

######  Proposition 3.2.

Let $S_{\beta}(t)$ and $P_{\beta}$ be the resolvent and integral resolvent
respectively generated by the operator $B$. Let $S_{\beta}^{\rho}(t)$ and
$P_{\beta}^{\rho}(t)$ be the resolvent and integral resolvent respectively
generated by the operator $-\rho I$. Then, we have the inequalities

(3.9) $S_{\beta}(t)\varphi(n)\leq S_{\beta}^{\rho}(t)\varphi(n)$

and

(3.10) $P_{\beta}(t)\varphi(n)\leq P_{\beta}^{\rho}(t)\varphi(n).$

###### Proof.

We use the representation (3.7) (3.8) and the fact that $B$ is the
infinitesimal generator of a Markovian semigroup $e^{tB}$

$\displaystyle S_{\beta}(t)\varphi(n)$
$\displaystyle=\int_{0}^{t}t^{-\beta}\Phi_{\beta}(\tau t^{-\beta})e^{\tau
B_{\rho}}\varphi(n)d\tau$
$\displaystyle=\int_{0}^{t}t^{-\beta}\Phi_{\beta}(\tau t^{-\beta})e^{\tau
B}e^{-t\rho}\varphi(n)d\tau$
$\displaystyle\leq\int_{0}^{t}t^{-\beta}\Phi_{\beta}(\tau
t^{-\beta})e^{-t\rho}\varphi(n)d\tau$
$\displaystyle=S_{\beta}^{\rho}(t)\varphi(n),$

Similarly for the integral resolvent, we consider any
$\varphi\in\mathcal{L}_{\gamma}$, we get

$\displaystyle P_{\beta}(t)\varphi(n)$
$\displaystyle=\int_{0}^{t}\beta\frac{\tau}{t^{1+\beta}}\Phi_{\beta}(\tau
t^{-\beta})e^{\tau B_{\rho}}\varphi(n)d\tau$
$\displaystyle=\int_{0}^{t}\beta\frac{\tau}{t^{1+\beta}}\Phi_{\beta}(\tau
t^{-\beta})e^{\tau B}e^{-t\rho}\varphi(n)d\tau$
$\displaystyle\leq\int_{0}^{t}\beta\frac{\tau}{t^{1+\beta}}\Phi_{\beta}(\tau
t^{-\beta})e^{-t\rho}\varphi(n)d\tau$
$\displaystyle=P_{\beta}^{\rho}(t)\varphi(n).$

∎

Next, we present the main result of this paper.

######  Theorem 3.3.

Let $\gamma>0$. Suppose the hypotheses of Lemma 3.1 holds and $f(x,\gamma)\leq
0$, for all $x\in\mathbb{R}$. Furthermore, suppose that $\partial_{s}f(x,0)$
is a non-increasing function of $x$. Let the initial value
$\varphi\in\mathcal{L}_{\gamma}$. Then, there exists a unique mild solution
$u\in C(\mathbb{R}_{0}^{+},\mathcal{L}_{\gamma})$ to the equation (3.6) and
hence a unique solution for nonlinear system (3.1).

###### Proof.

We used a iterative method. We define the following operator

(3.11)
$\mathcal{K}_{\beta}u(n,t):=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,u(n,s))\
ds,$

where $F(x,y)=\rho y+f(x,y)$ and inequality (3.3) holds. Note that
$F(x,0)=f(x,0)=0$. Furthermore, its follows from Dominated convergence
Lebesgue Theorem that the operator $\mathcal{K}_{\beta}$ is continuous over
$C(\mathbb{R}^{+},\mathcal{L}_{\gamma})$.

We construct the sequence

$v_{k}(n,t)=\begin{cases}\mathcal{K}_{\beta}v_{k-1}(n,t),&\ \mbox{if}\ k\geq
1,\\\ 0,&\ \mbox{if}\ k=0.\end{cases}$

We computing $v_{1}(n,t)$, obtaining

$v_{1}(n,t)=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,0)\
ds=S_{\beta}(t)\varphi(n)\geq 0,$

since $S_{\beta}(t)$ is a positive operator. Using the fact that $F(x,y)$ is
no negative for all $x\in\mathbb{R}$ and $0\leq y\leq\gamma$, we get

$\displaystyle v_{1}(n,t)$ $\displaystyle=S_{\beta}(t)\varphi(n)$
$\displaystyle\leq
S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,v_{1}(n,t))\ ds$
$\displaystyle=v_{2}(n,t).$

We remark that $F(x,y)$ is no decreasing in the second variable by Lemma
(3.1). Let $k\in\mathbb{N}$. We shall prove that $v_{k}(n,t)\leq
v_{k+1}(n,t)$. Indeed, we have

$\displaystyle v_{k}(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,v_{k}(n,t))\
ds$ $\displaystyle\leq
S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,v_{k+1}(n,t))\ ds$
$\displaystyle=v_{k+1}(n,t).$

Therefore, we get the inequalities $v_{0}\leq v_{1}\leq\ldots\leq v_{k}\leq
v_{k+1}\leq\ldots$.

On the other hand, we consider the sequence

$w_{k}(n,t)=\begin{cases}\mathcal{K}_{\beta}w_{k-1}(n,t),&\ \mbox{if}\ k\geq
1,\\\ \gamma,&\ \mbox{if}\ k=0.\end{cases}$

Then, using that $F(x,y)$ is no decreasing in the second variable, we get

$\displaystyle w_{1}(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,w_{0}(n,t))\
ds$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,\gamma)\
ds.$

Note the fact that $0\leq F(x,y)\leq\rho\gamma$ since $f(x,\gamma)\leq 0$.
Using this fact and the inequalities (3.9) and (3.10) we obtain

$\displaystyle w_{1}(n,t)$ $\displaystyle\leq
S_{\beta}^{\rho}(t)\varphi(n)+\int_{0}^{t}P_{\beta}^{\rho}(t-s)F(n-cs,\gamma)\
ds$ $\displaystyle\leq
S_{\beta}^{\rho}(t)\varphi(n)+\rho\gamma\int_{0}^{t}P_{\beta}^{\rho}(t-s)\
ds.$

We recall that $0\leq\varphi(n)\leq\gamma$ and the fact that
$S_{\beta}^{\rho}(0)=1$. Furthermore, we recall that

$\frac{d}{dt}S_{\beta}^{\rho}(t)=-\rho P_{\beta}^{\rho}(t).$

Therefore, integrating the above equality and computing, we have

$\displaystyle w_{1}(n,t)$ $\displaystyle\leq
S_{\beta}^{\rho}(t)\varphi(n)+\rho\gamma\int_{0}^{t}P_{\beta}^{\rho}(t-s)\ ds$
$\displaystyle\leq S_{\beta}^{\rho}(t)\gamma-S_{\beta}^{\rho}(t)\gamma+\gamma$
$\displaystyle=\gamma.$

We use above inequality and Lemma (3.1) for to obtain

$F(n-cs,w_{1}(n,t))\leq F(n-cs,\gamma).$

Hence, we obtain

$\displaystyle w_{1}(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,\gamma)\
ds$ $\displaystyle\leq
S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,w_{1})\ ds$
$\displaystyle=w_{2}(n,t).$

Let $k\in\mathbb{N}$. We shall prove that $w_{k+1}(n,t)\leq w_{k}(n,t)$.
Indeed, we have

$\displaystyle w_{k+1}(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,w_{k}(n,t))\
ds$ $\displaystyle\leq
S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,w_{k-1})\ ds$
$\displaystyle=w_{k}(n,t).$

Therefore, we get

$\gamma\geq w_{1}(n,t)\geq w_{2}(n,t)\geq\ldots\geq w_{k}(n,t)\geq
w_{k+1}(n,t)\geq\ldots.$

Assertion: For any $j,k\in\mathbb{N}_{0}$, we have

$0\leq v_{1}(n,t)\leq\ldots\leq v_{k}(n,t)\leq\ldots\leq
w_{j}(n,t)\leq\ldots\leq\omega_{1}(n,t)\leq\gamma.$

We suppose that there exists $j,k$ such that $w_{j}(n,t)\leq v_{k}(n,t)$.
Then, we can find the first entire number $m\in\mathbb{N}$ such that
$w_{m}(n,t)>v_{m}(n,t)$. Using lemma (3.1) we get

$\displaystyle 0<(v_{m}-w_{m})$
$\displaystyle=\int_{0}^{t}P_{\beta}(t-s)[F(n-cs,v_{m-1})-F(n-cs,w_{m-1})]\leq
0,$

which is a contradiction.

Therefore, it is clear that $\\{v_{k}\\}_{k\in\mathbb{N}}$ and
$\\{w_{k}\\}_{k\in\mathbb{N}}$ are convergent sequences. Let us define
$\displaystyle v:=\lim_{k\to\infty}v_{k}$ and $\displaystyle
w:=\lim_{k\to\infty}w_{k}$. It is clear from assertion that $v(n,t)\leq
w(n,t)$. Moreover, we have that $v,w\in
C(\mathbb{R}^{+},\mathcal{L}_{\gamma})$. Since the operator $K_{\beta}$ is
continuous over $C(\mathbb{R}^{+},\mathcal{L}_{\gamma})$, we have that $v$ and
$w$ are solution for (3.6).

We shall prove uniqueness of solution. Indeed, the application the Mean Value
Theorem yields

$\displaystyle 0\leq(w-v)(n,t)$
$\displaystyle=\int_{0}^{t}P_{\beta}(t-s)[F(n-cs,w(n,t))-F(n-cs,v(n,t))]$
$\displaystyle=\int_{0}^{t}P_{\beta}(t-s)\partial_{s}F(n-cs,\theta)[(w-v)(n,t)],$

where $v(n,t)\leq\theta\leq w(n,t)$. Since
$\partial_{s}F(n-cs,\theta)=\rho+\partial_{s}f(n-cs,\theta)$, the Lemma (3.1)
yields $\partial_{s}f(n-cs,\theta)$ is non-increasing and hence
$\partial_{s}f(n-cs,\gamma)\leq\partial_{s}f(n-cs,\theta)\leq\partial_{s}f(n-cs,0)$.
Consequently, we have
$\leq\partial_{s}F(n-cs,\theta)\leq\partial_{s}F(n-cs,0)$.

Observe that, for $y$ fixed, $\partial_{s}F(\cdot,y)$ is no decreasing. Then,
from above inequality, we obtain

$\displaystyle 0\leq(w-v)(n,t)$
$\displaystyle=\int_{0}^{t}P_{\beta}(t-s)[F(n-cs,w(n,t))-F(n-cs,v(n,t))]$
$\displaystyle=\partial_{s}F(n,0)\int_{0}^{t}P_{\beta}(t-s)[(w-v)(n,t)].$

Therefore, we get

$\displaystyle 0$ $\displaystyle\leq\sup_{n\in\mathbb{Z}_{+}}[w(n,t)-v(n,t)]$
$\displaystyle\leq\partial_{s}F(n,0)\int_{0}^{t}P_{\beta}(t-s)\sup_{n\in\mathbb{Z}_{+}}[(w-v)(n,t)]$
$\displaystyle\leq\partial_{s}F(n,0)\int_{0}^{t}P_{\beta}^{\rho}(t-s)\sup_{n\in\mathbb{Z}_{+}}[(w-v)(n,t)],$

where we used that (3.10). Finally, from Gronwall inequality (see e.g [42])
for integral inequalities, we obtain $v\equiv w$, getting the desired result.
∎

Next, we present two consequences of the Theorem (3.3) concerning to a
comparison principle between solution with respective initial values.

######  Corollary 3.4.

Assume the hypotheses of Theorem (3.3). Let $u_{1}(n,t)$ and $u_{2}(n,t)$ the
solutions for equation (3.6) with initial condition $\varphi$ and $\psi$
respectively. If $\varphi(n)\leq\psi(n)$, for all $n\in\mathbb{Z}$, then
$u_{1}(n,t)\leq u_{2}(n,t)$, for all $n\in\mathbb{Z}$ and $t\geq 0$.

###### Proof.

Analogously to proof of the Theorem (3.3), we consider the sequence
$\\{u_{1,k}\\}_{k\in\mathbb{N}}$ and $\\{u_{2,k}\\}_{k\in\mathbb{N}}$, defined
according to iterative definition (3). By uniqueness of solution for the
system (3.6), we have $u_{1,k}\to u_{1}$ and $u_{2,k}\to u_{2}$ when
$n\to\infty$. Since $u_{1,0}=u_{2,0}=0$ and $F(x,0)=0$, we have

$\displaystyle u_{1,1}(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,u_{1,0}(n,t))\
ds$ $\displaystyle=S_{\beta}(t)\varphi(n)$

and

$\displaystyle u_{2,1}(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,u_{2,0}(n,t))\
ds$ $\displaystyle=S_{\beta}(t)\psi(n).$

Since $S_{\beta}(t)$ is a positive operator, we have $u_{1,1}(n,t)\leq
u_{2,1}(n,t)$, for all $n\in\mathbb{Z}$ and $t\geq 0$.

Using Lemma (3.1), we get

$\displaystyle u_{1,2}(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,u_{1,1}(n,t))\
ds$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,u_{2,1}(n,t))\
ds$ $\displaystyle=u_{2,2}(n,t).$

By induction principle, we obtain $u_{1,k}(n,t)\leq u_{2,k}(n,t)$ for all
$k\in\mathbb{N}$ and therefore, $u_{1}(n,t)\leq u_{2}(n,t)$. The proof is
finished. ∎

######  Corollary 3.5.

Assume the hypotheses of Theorem (3.3). Let $v,w\in
C(\mathbb{R}_{0}^{+},\mathcal{L}_{\gamma})$ such that
$\mathcal{K}_{\beta}v(n,t)\geq v(n,t)$ and $\mathcal{K}_{\beta}w(n,t)\leq
w(n,t)$, for all $n\in\mathbb{Z}$ and $t\in\mathbb{R}_{0}^{+}$. Suppose that
$v(n,0)\leq w(n,0)$ for all $n\in\mathbb{Z}$. Then the inequality $v(n,t)\leq
w(n,t)$ for all $n\in\mathbb{Z}$ and $t\in\mathbb{R}_{0}^{+}$ holds.

###### Proof.

We use the Lemma (3.1) in order to obtain

$\displaystyle\mathcal{K}_{\beta}^{2}v(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,\mathcal{K}_{\beta}v(n,t))\
ds$ $\displaystyle\geq
S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,v(n,t))\ ds$
$\displaystyle=\mathcal{K}_{\beta}v(n,t).$

Applying recursively the above process, we get

$v(n,t)\leq\mathcal{K}_{\beta}v(n,t)\leq\mathcal{K}_{\beta}^{2}v(n,t)\leq\ldots\leq\lim_{k\to\infty}\mathcal{K}_{\beta}^{k}v(n,t):=v^{*}(n,t).$

Analogously, from Lemma (3.1), we obtain

$\displaystyle\mathcal{K}_{\beta}^{2}w(n,t)$
$\displaystyle=S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,\mathcal{K}_{\beta}w(n,t))\
ds$ $\displaystyle\leq
S_{\beta}(t)\varphi(n)+\int_{0}^{t}P_{\beta}(t-s)F(n-cs,w(n,t))\ ds$
$\displaystyle=\mathcal{K}_{\beta}w(n,t).$

and proceeding inductively, we get

$w(n,t)\geq\mathcal{K}_{\beta}w(n,t)\geq\mathcal{K}_{\beta}^{2}w(n,t)\geq\ldots\geq\lim_{k\to\infty}\mathcal{K}_{\beta}^{k}w(n,t):=w^{*}(n,t).$

Note that $v^{*}$ and $w^{*}$ there exists, since the operator
$\mathcal{K}_{\beta}$ is continuous. Furthermore, we remark that $v^{*}$ and
$w^{*}$ are solution for the equation (3.6) with initial values $v(n,0)$ and
$w(n,0)$ respectively. Then, the Corollary (3.4) yields

$v(n,t)\leq v^{*}(n,t)\leq w^{*}(n,t)\leq w(n,t),$

for all $n\in\mathbb{Z}$ and $t\geq 0$. The proof is finished. ∎

##  4\. Applications and Examples

In this section we present convolution operators and nonlinear functions
satisfying the hypotheses of Theorem (3.3). We conclude the paper with a
historical review about some cases and its applications in diverse models.

###  4.1. Convolution Operators

We present a sufficient condition in order to guarantee the Markov properties
for a $C_{0}$-semigroup generated by a discrete convolution operator.

######  Theorem 4.1.

Let $A\in\mathcal{B}(\ell^{p}(\mathbb{Z}))$ be an operator by discrete
convolution, namely $Af(n)=(a\ast f)(n)$, which is the infinitesimal generator
of Markovian semigroup $e^{tA}$ and $\alpha\in(0,1)$. Then, the semigroup
generated for the fractional powers $-(-A)^{\alpha}$, namely
$e^{-t(-A)^{\alpha}}$ is a Markovian semigroup.

###### Proof.

The key of the proof is the Lévy subordination principle and its properties.
Let $\phi\geq 0$ be a sequence in $\ell^{p}(\mathbb{Z})$, $1\leq p\leq\infty$.
Since $\alpha\in(0,1)$, from remark (5.2) property i) we have

$e^{-t(-A)^{\alpha}}\phi(n)=\int_{0}^{\infty}e^{\lambda
A}\phi(n)f_{t,\alpha}(\lambda)d\lambda,\quad t>0.$

The positivity of Lévy process (see remark (5.2), property ii)) and $e^{tA}$
(Markovian property) yield

$e^{-t(-A)^{\alpha}}\phi(n)\geq 0.$

Let us choose $\phi(n)=\delta_{0}(n)$. A straightforward computation yield

$0\leq\int_{0}^{\infty}e^{\lambda
A}\delta_{0}(n)f_{t,\alpha}(\lambda)d\lambda=e^{-t(-A)^{\alpha}}\delta_{0}(n)=e^{-t(-a)^{\alpha}}.$

Let $\textbf{1}(n)$ be the sequence constantly equals to $1$ considered in
(2.3). We use again the integral representation of Lévy subordination
principle, and (2.3) in order to obtain

$e^{-t(-A)^{\alpha}}\textbf{1}(n)=\int_{0}^{\infty}e^{\lambda
A}\textbf{1}(n)f_{t,\alpha}(\lambda)d\lambda=\textbf{1}(n)\int_{0}^{\infty}f_{t,\alpha}(\lambda)d\lambda=\textbf{1}(n).$

On the other hand, since $e^{-t(-A)^{\alpha}}$ is defined by discrete
convolution, we have

$e^{-t(-A)^{\alpha}}\textbf{1}(n)=\sum_{m\in\mathbb{Z}}e^{-t(-a)^{\alpha}}(m)\textbf{1}(n-m)=\textbf{1}(n)\sum_{m\in\mathbb{Z}}e^{-t(-a)^{\alpha}}(m).$

Therefore, from above equalities, we have
$\displaystyle\|e^{-(-a)^{\alpha}}\|_{1}=\sum_{n\in\mathbb{Z}}e^{-t(-a)^{\alpha}}(n)=1$.
From (2.2), we obtain $\|e^{-(-A)^{\alpha}}\|_{1}=1$. The proof is finished.

∎

######  Corollary 4.2.

Let $\alpha\in(0,1)$. The semigroups $e^{-t(-\Delta_{d})^{\alpha}}$,
$e^{-t(\Delta/-\nabla)^{\alpha}}$ and $e^{-t(-\Delta_{dd})^{\alpha}}$ are
Markovian.

###### Proof.

A straightforward computation yields that $e^{t\Delta_{d}}$ and
$e^{t(-\Delta/\nabla)}$ are Markovian semigroups. Indeed, from generatrix
formula for $I_{\nu}$ Bessel function, we have

$\|e^{t\Delta_{d}}\|_{1}=e^{-2t}\sum_{n\in\mathbb{Z}}I_{n}(2t)=1$

and

$\|e^{t\Delta_{dd}}\|_{1}=e^{-2t}\sum_{n\in\mathbb{Z}}I_{\frac{n}{2}}\chi_{2\mathbb{Z}}(2t)=1.$

Furthermore

$\|e^{t(-\Delta/\nabla)}\|_{1}=e^{-t}\sum_{n\in\mathbb{Z}}\frac{t^{(\pm
n)}}{(\pm n)!}\chi_{\mathbb{Z}_{\pm}}(n)=1.$

The positive condition is immediate from positivity of $I_{n}(2t),e^{-2t}$ and
$t^{n}$, since $t\geq 0$.

Therefore, from the above proposition, we obtain the desired result. ∎

######  Remark 4.3.

The above result extends the result obtained by Lizama and Roncal in [23],
Theorem 1.3, v) in the study of the discrete fractional Laplacian and the
properties of the generated semigroup.

###  4.2. Nonlinear Functions

Next, we present nonlinear functions which satisfies the hypotheses of Theorem
(3.3).

  1. (1)

Let $a\in\mathbb{R}^{+}$ a fixed constant and define

(4.1) $f(s)=s(a^{2}-s^{2}).$

It is clear that $f$ is a measurable function, the restriction $f$ belongs to
$C^{1}(0,a)$ and $f(0)=f(a)=0$. Its partial derivatives
$\partial_{s}f(s)=a^{2}-3s^{2}$ and $\partial_{s}^{2}f(s)=-6s<0$ gives us the
concavity of $f$ on $[0,a]$. Let us choose $\rho=3a^{2}$. With this choice, we
ensure that the inequality (3.3) holds. Furthermore, its clear that
$\partial_{s}f(x,0)$ is a no decreasing function, and
$f(x,\gamma)=\gamma(r(x)-\gamma^{p})=\gamma^{\frac{1}{p}}(r(x)-\gamma)\leq 0$.
Therefore, the function $f$ defined in (4.1) satisfies all the hypotheses of
the Theorem (3.3).

  2. (2)

Let $r:\mathbb{R}\to\mathbb{R}$ be continuous, no decreasing, bounded and
piecewise continuously differentiable satisfying $0<r(\infty)<\infty$. In
fact, we define $r(\infty)=\gamma$. We define the nonlinear function
$f:\mathbb{R}\times\mathbb{R}^{+}\to\mathbb{R}$ as follows

(4.2) $f(x,s)=s(r(x)-s).$

It is clear that $f$ is a measurable function, the restriction $f(x,\cdot)$
belongs to $C^{1}(0,\gamma)$ and $f(x,0)=0$, for all $x\in\mathbb{R}$. Its
partial derivatives $\partial_{s}f(x,s)=r(x)-2s$ and
$\partial_{s}^{2}f(x,s)=-2<0$ gives us the concavity of $f$ on $[0,\gamma]$.
Let us choose $\rho=3\gamma$. With this choice, we ensure that the inequality
(3.3) holds. Furthermore, its clear that $\partial_{s}f(x,0)$ is a no
decreasing function, and $f(x,\gamma)=\gamma(r(x)-\gamma)\leq 0$. Therefore,
the function $f$ defined in (4.2) satisfies all the hypotheses of the Theorem
(3.3).

  3. (3)

We can generalize the above function. Let $r:\mathbb{R}\to\mathbb{R}$ be
continuous, no decreasing, bounded and piecewise continuously differentiable
satisfying $0<r(\infty)<\infty$. In fact, we define
$r(\infty)=\gamma^{\frac{1}{p}}$. Let $p\geq 1$. We define the nonlinear
function $f:\mathbb{R}\times\mathbb{R}^{+}\to\mathbb{R}$ as follows

(4.3) $f(x,s)=s(r(x)-s^{p}).$

It is clear that $f$ is a measurable function, the restriction $f(x,\cdot)$
belongs to $C^{1}(0,\gamma)$ and $f(x,0)=0$, for all $x\in\mathbb{R}$. Its
partial derivatives $\partial_{s}f(x,s)=r(x)-(p+1)s^{p}$ and
$\partial_{s}^{2}f(x,s)=-p(p+1)s^{p-1}<0$ gives us the concavity of $f$ on
$[0,\gamma]$. Let us choose $\rho=(p+2)\gamma^{p}$. With this choice, we
ensure that the inequality (3.3) holds. Furthermore, its clear that
$\partial_{s}f(x,0)$ is a no decreasing function, and
$f(x,\gamma)=\gamma(r(x)-\gamma^{p})=\gamma^{\frac{1}{p}}(r(x)-\gamma)\leq 0$.
Therefore, the function $f$ defined in (4.3) satisfies all the hypotheses of
the Theorem (3.3).

######  Remark 4.4.

The model (3.1), considering the fractional Laplacian and the nonlinear
function (4.1) is the well-known _discrete fractional generalized Fisher-KPP
equation_. This includes the discrete generalized Fisher-KPP equation with
delay since $f(n-ct,s)=s(r(n-ct)-s)$ with $n\in\mathbb{Z}$ and some constant
$c>0$. The classical Fisher-KPP equation corresponds to the case where
$r(x)\equiv\gamma$.

When we considerate $r\equiv 1$ and $p=2$, we recover the _Newell-Whitehead-
Segel equation_ , which describes the so called Rayleigh-Benard convection.
The Newell-Whitehead-Segel equation is a well-known universal equation to
govern evolution of nearly one-dimensional nonlinear patterns produced by a
finite-wavelength instability in isotropic two-dimensional media. More
references for these applications, see [26] and [31].

##  5\. Appendix

In this appendix we give the basic concepts of the special functions involving
in the development of this study.

###  5.1. Wright function $\Phi_{\gamma}$

The Wright type function with one parameter from [13, Formula (28)] (see also
[28, 30, 41]) is given by

(5.1) $\Phi_{\alpha}(z):=\sum_{n=0}^{\infty}\frac{(-z)^{n}}{n!\Gamma(-\alpha
n+1-\alpha)}=\frac{1}{2\pi
i}\int_{\gamma}\mu^{\alpha-1}e^{\mu-z\mu^{\alpha}}d\mu,\quad 0<\alpha<1,$

where $\gamma$ is a contour which starts and ends at $-\infty$ and encircles
the origin once counterclockwise. This special case has sometimes been called
the Mainardi function.

An interesting fact is the following relationship between the Wright and the
Mittag-Leffler functions.

######  Remark 5.1.

Let $z\in\mathbb{C}$, $t>0$ and $0<\alpha,\gamma<1$. Then the following
properties hold:

  * (i)

$\displaystyle E_{\gamma,1}(z)=\int_{0}^{\infty}\Phi_{\gamma}(t)e^{zt}dt$.

  * (ii)

$\displaystyle\Phi_{\alpha}(t)\geq 0$.

  * (iii)

$\displaystyle\int_{0}^{\infty}\Phi_{\alpha}(t)dt=1$.

It follows from (ii) and (iii) that $\Phi_{\alpha}$ is a probability density
function on $\mathbb{R}_{+}$. Actually, the Wright function has been used for
models in stochastic processes [13, 14].

###  5.2. Lévy process $f_{t,\alpha}$

We present the following function, called stable Lévy process, defined for
$0<\alpha<1$ by

(5.2) $f_{t,\alpha}(\lambda)=\begin{cases}\displaystyle\frac{1}{2\pi
i}\int_{\sigma-i\infty}^{\sigma+i\infty}e^{z\lambda-
tz^{\alpha}}dz,\;\;\;\;\sigma>0,\quad t>0,&\lambda\geq 0,\\\
0\quad&\lambda<0,\end{cases}$

where the branch of $z^{\alpha}$ is taken so that $\mbox{Re}(z^{\alpha})>0$
for $\mbox{Re}(z)>0.$ This branch is single-valued in the $z$-plane cut along
the negative real axis. These functions were introduced by S. Bochner [7] in
the study of certain stochastic processes. K. Yosida [43] used them
systematically in the study of semigroups generated by fractional powers of
uniformly bounded semigroups of linear operators. The Lévy functions are the
density functions associated with the stable Lévy processes in the rotational
invariant case, and are related to the fractional Brownian motion.

######  Remark 5.2.

The following properties hold:

  * (i)

$\displaystyle\int_{0}^{\infty}e^{-\lambda
a}f_{t,\alpha}(\lambda)d\lambda=e^{-ta^{\alpha}},\quad t>0,\quad a>0,\quad
0<\alpha<1.$

  * (ii)

$\displaystyle f_{t,\alpha}(\lambda)\geq 0,\quad\lambda>0,\,\,t>0,\quad
0<\alpha<1.$

  * (iii)

$\displaystyle\int_{0}^{\infty}f_{t,\alpha}(\lambda)d\lambda=1,\quad t>0,\quad
0<\alpha<1.$

For a proof of (i)-(iii), see [43, p.260-262].

###  5.3. Bessel function $I_{\nu}$

For $\nu\in\mathbb{R}$, the _Modified Bessel functions of the first kind_ is
defined by

(5.3)
$I_{\nu}(x)=\sum_{n=0}^{\infty}\frac{1}{\Gamma(n+\nu+1)n!}\left(\frac{x}{2}\right)^{2n+\nu}.$

It is direct from the definition (5.3) that $I_{n}(x)\geq 0,\quad
n\in\mathbb{Z},\,\,x\geq 0$, Following propierties can be founded in [15,
Formula 8.511], [4].

  1. (1)

$\displaystyle\sum_{n\in\mathbb{Z}}I_{n}(x)z^{n}=\displaystyle
e^{\frac{x}{2}(z+\frac{1}{z})},\quad z\in\mathbb{C}\setminus\\{0\\}$.

  2. (2)

$I_{-n}(x)=I_{n}(x)=(-1)^{n}I_{n}(-x)$.

  3. (3)

$\displaystyle I_{n}(x+y)=\sum_{k\in\mathbb{Z}}I_{n-k}(x)I_{k}(y)$.

## References

  * [1] L. Abadías, J. González-Camus, P.J. Miana and J.C. Pozo. _Large time behaviour for the heat equation on $\mathbb{Z}$; moments and decay rates._ J. Math. Anal. Appl. 500 (2021) Issue 2 125137.
  * [2] L. Abadias, M. de León-Contreras and J.L. Torrea. _Non-local fractional derivatives. Discrete and continuous._ J. Math. Anal. Appl. 449(1) (2017). 734–755.
  * [3] L. Abadias and P. J. Miana. _A subordination principle on Wright functions and regularized resolvent families._ Journal of Function Spaces. Article ID 158145, 2015.
  * [4] G.E. Andrews, R. Askey and R. Roy._Special functions. Encyclopedia of mathematics and its applications._ Cambrigde University Press, Vol. 71, 2006.
  * [5] H. Bateman. _Some simple differential difference equations and the related functions._ Bull. Amer. Math. Soc. 49 (1943), 494–512.
  * [6] E. Bazhlekova. _Fractional evolution equations in Banach spaces._ Ph.D Thesis, Eindhoven University of Technology, 2001.
  * [7] S. Bochner. _Diffusion equation and stochastic processes._ Proc. Nat. Acad. Sci. U. S. A. 35 (1949), 368–370.
  * [8] Ó. Ciaurri, C. Lizama, L. Roncal and J. L. Varona. _On a connection between the discrete fractional Laplacian and superdiffusion._ Appl. Math. Letters. 49 (2015), 119–125.
  * [9] Ó. Ciaurri, T. A. Gillespie, L. Roncal, J. L. Torrea and J. L. Varona. _Harmonic Analysis associated with a discrete Laplacian._ J. Anal. Math. 132 (2017), 109–131.
  * [10] K.Engel and R. Nagel. _One-Parameters semigroups for linear evolution equations._ Springer-Verlag, Berlin. 1980.
  * [11] J. González-Camus, V. Keyantuo, C. Lizama and M. Warma. _Fundamental solutions for discrete dynamical systems involving the fractional Laplacian_. Math. Meth. Appl. Sci, (2019), 4688–4711.
  * [12] J. González-Camus, C. Lizama and P.J. Miana. _Fundamental solutions for semidiscrete evolution equations via Banach algebras_. Adv. Differ. Equ, 35 (2021).
  * [13] R. Gorenflo, Y. Luchko and F. Mainardi. _Analytical properties and applications of the Wright function._ Fract. Calc. Appl. Anal. 2 (1999), 383–414.
  * [14] R. Gorenflo and F. Mainardi. _On Mittag-Leffler-type functions in fractional evolution processes._ J. Comp. Appl. Math. 118 (2000), 283–299.
  * [15] I.S. Gradshteyn and I.M. Ryzhik. _Table of Integrals, Series and Products._ Seventh edition. Elsevier, Academic Press. Burlington, San Diego, London. 2007.
  * [16] C. Hu and B. Li. _Spatial dynamics for lattice differential equations with a shifting habitat._ J. Differential Equations. (2015) 1967–1989.
  * [17] J.P. Keener. _On the formation of circulation patterns of excitation in anisotropic excitable media._ J. Math. Biol. 26 (1988), 41–56.
  * [18] V. Keyantuo, C. Lizama and M. Warma. _Existence, regularity and representation of solutions of time fractional diffusion equations._ Adv. Differ. Equ. 21 (9-10) (2016), 837–886.
  * [19] V. Keyantuo, C. Lizama and M. Warma. _Spectral criteria for solvability of boundary value problems and positivity of solutions of time-fractional differential equations_. Abstract and Applied Analysis, vol. 2013, Article ID 614328, (2013), 11 pages.
  * [20] V. Keyantuo, C. Lizama and M. Warma. _Lattice dynamical systems associated with a fractional Laplacian._ Numer. Funct. Anal. Optim. 40 (2019), no. 11, 1315–1343.
  * [21] R. Larsen. _Banach Algebras: An Introduction._ Marcel Dekker, New York, 1973.
  * [22] C. Lizama and M. Murillo. _On a connection between the N-dimensional fractional Laplacian and 1-D operators on lattices._ J. Math. Anal. Appl. 511 (1) (2022), 126051.
  * [23] C. Lizama and L. Roncal. _Hölder-Lebesgue regularity and almost periodicity for semidiscrete equations with a fractional Laplacian._ Discrete Contin. Dyn. Syst, Series A, 38 (3)(2018), 1365–1403.
  * [24] F. Mainardi, Y. Luchko and G. Pagnini. _The fundamental solution of the space-time fractional diffusion equation._ Fract. Calc. Appl. Anal. 4 (2) (2001), 153–192.
  * [25] M. Ortigueira. _Fractional central differences and derivatives._ J. Vib. Control 14 (2008), no. 9–10, 1255–1266.
  * [26] A.C. Newell and J. A. Whitehead. _Finite bandwith, finite amplitude convection._ J. Fluid Mech. 38 (1969), 279–303.
  * [27] A Pazy. _Semigroups of linear operators and applications to partial differential equations._ Springer-Verlag, New York. 1983.
  * [28] I. Podlubny. _Fractional Differential Equations._ Academic Press, San Diego, CA, 1999.
  * [29] J. Prüss. _Evolutionary Integral Equations and Applications_. Birkhäuser, Basel, 1993.
  * [30] G. Samko, A. Kilbas and O. Marichev. _Fractional Integrals and Derivatives: Theory and Applications._ Gordon and Breach, New York and London, 1993.
  * [31] L.A. Segel. _Distant sidewalls cause slow amplitude modulation of cellular convection._ J. Fluid Mech. 38 (1969), 203–224.
  * [32] A. Slavík. _Asymptotic behavior of solutions to the semidiscrete diffusion equation._ Appl. Math. Lett. 106 (2020), 106392.
  * [33] A. Slavík and P. Stehlík. _Explicit solutions to dynamic diffusion-type equations and their time integrals_. Appl. Math. Comput., 234 (2014), 486–505.
  * [34] A. Slavík and P. Stehlík. _Dynamic diffusion-type equations on discrete-space domains._ J. Math. Anal. Appl. Vol 427, Issue 1, (2015), 525–545.
  * [35] P. Stehlík and J. Volek. _Transport equation on semidiscrete domains and Poisson–Bernoulli processes._ J. Difference Equ. Appl., 19 (3) (2013), 439–456.
  * [36] V.E. Tarasov. _Lattice fractional calculus._ Appl Math. Comp. 257 (2015), 12–33.
  * [37] V.E. Tarasov. _Fractional-order difference equations for physical lattices and some applications._ J Math Phys. 56 (2015), 103506.
  * [38] V.E. Tarasov. _Exact discretization of fractional Laplacian_. Comp. Math. Appl. 73 (2017) 855–863.
  * [39] O. Vallée and M. Soares. _Airy functions and applications to physics._ Second edition. Imperial College Press, London, 2010.
  * [40] J. Volek. _Maximum and minimum principles for nonlinear transport equations on discrete-space domains._ Electron. J. Differential Equations, 2014 (78) (2014), 1–13.
  * [41] E.M. Wright. _The generalized Bessel function of order greater than one._ Quart. J. Math., Oxford Ser.11 (1940), 36–48.
  * [42] H. Ye, J. Gao and Y. Ding. _A generalized Gronwall inequality and its application to a fractional differential equation._ J. Math. Anal. Appl. 328 (2) (2007), 1075–-1081.
  * [43] K. Yosida. _Functional Analysis._ Springer Verlag, Berlin Heidelberg New York, 1980.
  * [44] B. Zinner, G. Harris and W. Hudson. _Traveling wave fronts for the discrete Fisher’s equation._ J. Differential Equations. 105 (1993), 46–62.