Bayesian Regularization on Function Spaces via Q-Exponential Process

Shiwei Lan Shuyi Li Michael O’Connor

School of Mathematical & Statistical Sciences, Arizona State University,
Tempe, AZ 85287

###### Abstract

Regularization is one of the most important topics in optimization, statistics
and machine learning. To get sparsity in estimating a parameter
$u\in{\mathbb{R}}^{d}$, an $\ell_{q}$ penalty term, $\|u\|_{q}$, is usually
added to the objective function. What is the probabilistic distribution
corresponding to such $\ell_{q}$ penalty? What is the correct stochastic
process corresponding to $\|u\|_{q}$ when we model functions $u\in L^{q}$?
This is important for statistically modeling large dimensional objects, e.g.
images, with penalty to preserve certainty properties, e.g. edges in the
image. In this work, we generalize the $q$-exponential distribution (with
density proportional to) $\exp{(-\frac{1}{2}|u|^{q})}$ to a stochastic process
named _$Q$ -exponential (Q-EP) process_ that corresponds to the $L_{q}$
regularization of functions. The key step is to specify consistent
multivariate $q$-exponential distributions by choosing from a large family of
elliptic contour distributions. The work is closely related to Besov process
which is usually defined by the expanded series. Q-EP can be regarded as a
definition of Besov process with explicit probabilistic formulation and direct
control on the correlation length. From the Bayesian perspective, Q-EP
provides a flexible prior on functions with sharper penalty ($q<2$) than the
commonly used Gaussian process (GP). We compare GP, Besov and Q-EP in modeling
time series and reconstructing images and demonstrate the advantage of the
proposed methodology.

##  1 INTRODUCTION

Regularization on function spaces is one of the fundamental questions in
statistics and machine learning. Large dimensional objects such as images can
be viewed as evaluations of proper functions. Statistical models for these
objects on function spaces demand regularization to induce sparsity, prevent
over-fitting, and produce meaningful reconstruction etc. Gaussian process [GP
22, 10] has been widely used as an $L_{2}$ penalty or a prior on the function
space. Despite of the flexibility, sometimes random candidate functions drawn
from GP are over-smooth for modeling certain objects such as images with sharp
edges. To address this issue, researchers have proposed a class of $L_{1}$
penalty based priors including Laplace random field [21, 17, 14] and Besov
process [16, 5, 11, 6]. They have been extensively applied in spatial modeling
[21], signal processing [14], imaging analysis [24, 17] and inverse problems
[16, 5]. Figure 1 demonstrates the random draws from GP, Besov and our
proposed $q$-exponential process respectively. The sample function by GP is
very smooth that contains large continuous blocks while samples from Besov and
$q$-exponential process manifest more pieces with “edges".

Given an orthonormal basis $\\{\phi_{\ell}\\}_{\ell=1}^{\infty}$ in
$L^{2}(D)$, with $D\subset{\mathbb{R}}^{d^{\star}}$ being the domain of
interest, a random function $u$ from the Besov process can be represented in
the following series [5]:

$u(x)=\sum_{\ell=1}^{\infty}\gamma_{\ell}u_{\ell}\phi_{\ell}(x),\quad
u_{\ell}\overset{iid}{\sim}\pi_{q}(\cdot)\propto\exp{(-\frac{1}{2}|\cdot|^{q})}$
(1)

where $q\geq 1$ and
$\gamma_{\ell}=\kappa^{-\frac{1}{q}}\ell^{-(\frac{s}{d^{\star}}+\frac{1}{2}-\frac{1}{q})}$
with (inverse) variance $\kappa>0$ and smoothness $s>0$. When $q=2$, this
reduces to GP but Besov is often used with $q=1$ to provide “edge-preserving"
function candidates suitable for image analysis. Though straightforward, such
series definition lacks a direct way to specify the correlation length as GP
does through the covariance function. What is more, once the basis
$\\{\phi_{\ell}\\}$ is chosen, there is no natural way to make prediction with
Besov process.

Figure 1: Sample functions from three priors: GP (left), Besov (middle) and
Q-EP (right).

We propose a novel stochastic process named _$q$ -exponential process (Q-EP)_
to address these issues. we start with the $q$-exponential distribution
$\pi_{q}(\cdot)$ and generalize it to a multivariate distribution (from a
family of elliptic contour distributions) that is consistent to
marginalization. Such consistency requires the joint distribution and the
marginalized one (by any subset of components) to have the same format of
density (See Section 3). We further generalize such multivariate
$q$-exponential distribution to the process Q-EP. The proposed Q-EP process is
related to the elliptic process in [1] which however focuses on a completely
different process named squeezebox process. Our work on Q-EP makes multi-fold
contributions to the functional regularization for statistical modeling:

  1. 1.

We propose a stochastic process Q-EP corresponding to the $L_{q}$
regularization on function spaces.

  2. 2.

We provide an explicit probabilistic formulation of Besov process through Q-EP
with direct ways to specify the correlation length and to make prediction.

  3. 3.

We prove that the elliptic slice sampler [ESS 18] can be applied to the
inference of Q-EP based models.

The rest of the paper is organized as follows. Section 2 introduces the
$q$-exponential distribution and its multivariate generalizations. We propose
the Q-EP with details in Section 3. In Section 4 we show that ESS can be
applied to models with Q-EP priors and demonstrate the advantage of Q-EP over
GP and Besov in time series modeling and image reconstruction. Finally we
discuss some future directions in Section 5.

##  2 THE $Q$-EXPONENTIAL DISTRIBUTION AND ITS MULTIVARIATE GENERALIZATIONS

Let us start with the $q$-exponential distribution for a scalar random
variable $u\in{\mathbb{R}}$. It is named in [5] and defined with the following
density not in an exact form (as a probability density normalized to 1):

$\pi_{q}(u)\propto\exp{(-\frac{1}{2}|u|^{q})}.$ (2)

This $q$-exponential distribution (2) is actually a special case of the
following _exponential power (EP)_ distribution $\mathrm{EP}(\mu,\sigma,q)$
with $\mu=0$, $\sigma=1$:

$p(u|\mu,\sigma,q)=\frac{q}{2^{1+1/q}\sigma\Gamma(1/q)}\exp\left\\{-\frac{1}{2}\left|\frac{u-\mu}{\sigma}\right|^{q}\right\\}$
(3)

Note the parameter $q>0$ in (3) controls the tail behavior of the
distribution: the smaller $q$ the heavier tail and vice verse. This
distribution also includes many commonly used ones such as the normal
distribution ${\mathcal{N}}(\mu,\sigma^{2})$ for $q=2$ and the Laplace
distribution $L(\mu,b)$ with $\sigma=2^{-1/q}b$ when $q=1$.

How can we generalize it to multivariate distribution and further to a
stochastic process? Gomez [9] provided one possibility of a multivariate EP
distribution, denoted as $\mathrm{EP}_{d}(\bm{\mu},{\bf C},q)$, with the
following density:

$\displaystyle p({\bf u}|\bm{\mu},{\bf C},q)=$
$\displaystyle\frac{q\Gamma(\frac{d}{2})}{2\Gamma(\frac{d}{q})}2^{-\frac{d}{q}}\pi^{-\frac{d}{2}}|{\bf
C}|^{-\frac{1}{2}}\cdot$ (4)
$\displaystyle\exp\left\\{-\frac{1}{2}\left[{({\bf
u}-\bm{\mu})}^{\mathsf{T}}{\bf C}^{-1}({\bf
u}-\bm{\mu})\right]^{\frac{q}{2}}\right\\}$

When $q=2$, it reduces to the familiar multivariate normal (MVN) distribution
${\mathcal{N}}_{d}(\bm{\mu},{\bf C})$, which is the foundation of Gaussian
process (GP), defined as a collection of random variables whose (any) finite
collection follows an MVN.

Unfortunately, the Gomez’s EP distribution $\mathrm{EP}_{d}(\bm{\mu},{\bf
C},q)$ cannot be generalized to a valid stochastic process because it does not
satisfy the marginalization consistency as MVN does (See Section 3 for more
details). It turns out we need to seek candidates in an even larger family
named _elliptic_ (contour) distribution $\mathrm{EC}_{d}(\bm{\mu},{\bf C},g)$:

######  Definition 1 (Elliptic (contour) distribution).

A multivariate elliptic distribution has the following density [12]

$p({\bf u})=k_{d}|{\bf C}|^{-\frac{1}{2}}g(r),\quad r({\bf u})={({\bf
u}-\bm{\mu})}^{\mathsf{T}}{\bf C}^{-1}({\bf u}-\bm{\mu})$ (5)

where $k_{d}>0$ is the normalizing constant and $g(\cdot)$, a one-dimensional
real-valued function independent of $d$ and $k_{d}$, is named _density
generating function_ [8].

Every random vector ${\bf u}\sim\mathrm{EC}_{d}(\bm{\mu},{\bf C},q)$ has a
stochastic representation [23, 3, 12]:

${\bf u}=\bm{\mu}+R{\bf L}S$ (6)

where $S\sim\mathrm{Unif}({\mathcal{S}}^{d+1})$ uniformly distributed on the
unit-sphere ${\mathcal{S}}^{d+1}$, ${\bf L}$ is the Cholesky factor of ${\bf
C}$ such that ${\bf C}={\bf L}{{\bf L}}^{\mathsf{T}}$, and $R^{2}=r({\bf
u})\sim
f(r)=\frac{\pi^{\frac{d}{2}}}{\Gamma(\frac{d}{2})}k_{d}r^{\frac{d}{2}-1}g(r)$.

The Gomez’s EP distribution $\mathrm{EP}_{d}(\bm{\mu},{\bf C},q)$ is a special
elliptic distribution $\mathrm{EC}_{d}(\bm{\mu},{\bf C},g)$ with
$g(r)=\exp\\{-\frac{1}{2}r^{\frac{q}{2}}\\}$. Not all elliptical distributions
can be used to create a valid process [1]. In the following, we will carefully
choose the density generator $g$ in $\mathrm{EC}_{d}(\bm{\mu},{\bf C},g)$ to
define a consistent multivariate $q$-exponential distribution ready to be
generalized to a process appropriately.

##  3 THE $Q$-EXPONENTIAL PROCESS

To generalize $\mathrm{EC}_{d}(\bm{\mu},{\bf C},g)$ to a valid stochastic
process, we need to choose proper $g$ such that the resulting distribution
satisfies two conditions of Kolmogorov’ extensition theorem [20]:

  1. 1.

Exchangeability The joint distribution needs to be invariant under finite
permutations, i.e. for any finite $d$ and permutation, $p({\bf
u}_{1:d})=p({\bf u}_{\pi(1:d)})$.

  2. 2.

Consistency The underlying distribution must be consistent, which means that
the marginal distribution of any collection of the random variables belongs to
the same distribution family as the original distribution, i.e. $p({\bf
u}_{1})=\int p({\bf u}_{1},{\bf u}_{2})d{\bf u}_{2}$.

As pointed out by Kano [13], the elliptic distribution
$\mathrm{EC}_{d}(\bm{\mu},{\bf C},g)$ in the format of Gomez’s EP distribution
(4) with $g(r)=\exp\\{-\frac{1}{2}r^{\frac{q}{2}}\\}$ does not satisfy the
consistency condition [also c.f. Proposition 5.1 of 9]. However, the following
theorem [13] suggests a different viable choice of $g$ to make a valid
generalization of $\mathrm{EC}_{d}(\bm{\mu},{\bf C},g)$ to a stochastic
process [1]:

######  Theorem 3.1.

An elliptic distribution is _consistent_ if and only if its density generator
function, $g(\cdot)$, has the following form

$g(r)=\int_{0}^{\infty}\left(\frac{s}{2\pi}\right)^{\frac{d}{2}}\exp\left\\{-\frac{rs}{2}\right\\}p(s)ds$
(7)

where $p(s)$ is a strictly positive mixing distribution independent of $d$ and
$p(s=0)=0$.

###  3.1 Consistent Multivariate $Q$-exponential Distribution

In the above theorem 3.1, if we choose $p(s)=\delta_{r^{\frac{q}{2}-1}}(s)$,
then we can have

$g(r)=r^{\left(\frac{q}{2}-1\right)\frac{d}{2}}\exp\left\\{-\frac{r^{\frac{q}{2}}}{2}\right\\}$
(8)

in (5) to define the following consistent _multivariate $q$-exponential
distribution_ $\mathrm{q}\\!-\\!\mathrm{ED}(\bm{\mu},{\bf C})$.

######  Definition 2.

A multivariate $q$-exponential distribution, denoted as
$\mathrm{q}\\!-\\!\mathrm{ED}(\bm{\mu},{\bf C})$, has the following density

$\displaystyle p({\bf u}|\bm{\mu},{\bf C},q)$
$\displaystyle=\frac{q}{2}(2\pi)^{-\frac{d}{2}}|{\bf
C}|^{-\frac{1}{2}}\boxed{r^{(\frac{q}{2}-1)\frac{d}{2}}}\exp\left\\{-\frac{r^{\frac{q}{2}}}{2}\right\\},$
(9) $\displaystyle r({\bf u})$ $\displaystyle={({\bf
u}-\bm{\mu})}^{\mathsf{T}}{\bf C}^{-1}({\bf u}-\bm{\mu})$

Regardless of the normalizing constant, our proposed multivariate
$q$-exponential distribution $\mathrm{q}\\!-\\!\mathrm{ED}(\bm{\mu},{\bf C})$
differs from the Gomez’s EP distribution $\mathrm{EP}_{d}(\bm{\mu},{\bf C},q)$
by a boxed term $r^{(\frac{q}{2}-1)\frac{d}{2}}$.

The following proposition determines the distribution of $R=\sqrt{r({\bf u})}$
thus gives a complete recipe for generating random vector ${\bf
u}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}(0,{\bf C})$ based on the stochastic
representation (6).

######  Proposition 3.1.

If ${\bf u}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}(0,{\bf C})$, then we have

$R^{q}=r^{\frac{q}{2}}\sim\Gamma\left(\alpha=\frac{d}{2},\beta=\frac{1}{2}\right)=\chi^{2}_{d}$
(10)

###### Proof.

With $g(r)$ chosen as in (8), the density of $r$ becomes

$f(r)\propto
r^{\frac{d}{2}-1}r^{\left(\frac{q}{2}-1\right)\frac{d}{2}}\exp\left\\{-\frac{r^{\frac{q}{2}}}{2}\right\\}=r^{\frac{q}{2}\cdot\frac{d}{2}-1}\exp\left\\{-\frac{r^{\frac{q}{2}}}{2}\right\\}$

A change of variable completes the proof. ∎

######  Corollary 3.1.

For each $k\in{\mathbb{N}}$, the $k$-th moment of $R$ is

$\mathrm{E}[R^{k}]=2^{\frac{k}{q}}\frac{\Gamma(\frac{d}{2}+\frac{k}{q})}{\Gamma(\frac{d}{2})}\overset{\cdot}{\sim}d^{\frac{k}{q}},\quad
as\quad d\to\infty$ (11)

###### Proof.

Since $v:=R^{q}\sim\Gamma\left(\alpha=\frac{d}{2},\beta=\frac{1}{2}\right)$ by
Proposition 3.1, we have:

$\begin{split}\mathrm{E}[R^{k}]&=\frac{1}{\Gamma(\frac{d}{2})}\left(\frac{1}{2}\right)^{\frac{d}{2}}\int_{0}^{\infty}v^{\frac{k}{q}+\frac{d}{2}-1}\exp\left\\{-\frac{1}{2}v\right\\}dv\\\
&=2^{\frac{k}{q}}\frac{\Gamma(\frac{d}{2}+\frac{k}{q})}{\Gamma(\frac{d}{2})}\overset{\cdot}{\sim}2^{\frac{k}{q}}\left(\frac{d}{2}\right)^{\frac{k}{q}}=d^{\frac{k}{q}}\end{split}$

where we use $\Gamma(x+\alpha)\overset{\cdot}{\sim}\Gamma(x)x^{\alpha}$ as
$x\to\infty$ with $x=\frac{d}{2}$ and $\alpha=\frac{k}{q}$ when $d\to\infty$.
∎

The following proposition [c.f. Theorem 2.6.4 in 8] tells the role of matrix
${\bf C}$ in characterizing the covariance between the components.

######  Proposition 3.2.

If ${\bf u}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}(\bm{\mu},{\bf C})$, then we
have

$\mathrm{E}[{\bf u}]=\bm{\mu},\;\mathrm{Cov}({\bf
u})=\frac{2^{\frac{2}{q}}\Gamma(\frac{d}{2}+\frac{2}{q})}{d\Gamma(\frac{d}{2})}{\bf
C}\overset{\cdot}{\sim}d^{\frac{2}{q}-1}{\bf C},\;as\;d\to\infty$ (12)

###### Proof.

By Theorem 2.6.4 in [8] for $\mathrm{q}\\!-\\!\mathrm{ED}_{d}(\bm{\mu},{\bf
C})=\mathrm{EC}_{d}(\bm{\mu},{\bf C},g)$ with $g$ chosen in (8), we know
$\mathrm{E}[{\bf u}]=\bm{\mu}$ and $\mathrm{Cov}({\bf
u})=(\mathrm{E}[R^{2}]/\textrm{rank}({\bf C})){\bf C}$. It follows by letting
$k=2$ in Corollary 3.1 and using the similar asymptotic analysis. ∎

To generalize ${\bf u}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}(0,{\bf C})$ to a
stochastic process, we want to scale it to ${\bf
u}^{*}=d^{\frac{1}{2}-\frac{1}{q}}{\bf u}$ so that we have a finite asymptotic
covariance. If ${\bf u}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}(0,{\bf C})$, then
we denote ${\bf u}^{*}\sim\mathrm{q}^{*}\\!-\\!\mathrm{ED}_{d}(0,{\bf C})$ as
_scaled_ $q$-exponential distribution. Now we are ready to define the _$q$
-exponential process (Q-EP)_ with the scaled $q$-exponential distribution.

######  Definition 3 (Q-EP).

A (centered) $q$-exponential process $u(x)$ with kernel ${\mathcal{C}}$,
$\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$, is a collection of random
variables such that any finite set, ${\bf u}=(u(x_{1}),\cdots u(x_{d}))$,
follows a scaled multivariate $q$-exponential distribution, i.e. ${\bf
u}\sim\mathrm{q}^{*}\\!-\\!\mathrm{ED}_{d}(0,{\bf C})$.

######  Remark 1.

When $d=1$, if we let $C=1$, then we have the density for $u$ as
$p(u)\propto|u|^{\frac{q}{2}-1}\exp\left\\{-\frac{1}{2}|u|^{q}\right\\}$,
differing from the original un-normalized density $\pi_{q}$ in (2) by a term
$|u|^{\frac{q}{2}-1}$. This is needed for the consistency of process
generalization. Numerically, it is dominated by the original
$\exp\left\\{-\frac{1}{2}|u|^{q}\right\\}$ and does not affect much of the
“edge-preserving" property of the Besov prior.

###  3.2 Connection and Contrast to the Besov Process

Both Besov and Q-EP are valid stochastic processes stemming from the
$q$-exponential distribution $\pi_{q}$. They are both designed to generalize
GP to have sharper regularization (through $q$) but Q-EP has advantages in
specifying correlation structure and making prediction.

It follows immediately that the covariance of the Besov process $u(\cdot)$ at
two points $x,x^{\prime}\in{\mathbb{R}}^{d^{\star}}$ is given

$\mathrm{Cov}(u(x),u(x^{\prime}))=\sum_{\ell=1}^{\infty}\gamma_{\ell}^{2}\phi_{\ell}(x)\otimes\phi_{\ell}(x^{\prime})$
(13)

Compared with (12), we have less control on the correlation length once the
orthonormal basis $\\{\phi_{\ell}\\}$ is chosen. On the other hand, Q-EP has
more freedom on the correlation structure through (12) with flexible choices
from a large class of kernels including powered exponential, Matérn, etc.
where we can directly specify the correlation length.

On the other hand, the following theorem states that Q-EP can also have series
expansion comparable to (1) for Besov.

######  Theorem 3.2 (Karhunen-Loéve).

If $u(x)\sim\mathrm{q}\\!-\\!\mathcal{EP}(0,{\mathcal{C}})$ and
${\mathcal{C}}$ is a trace operator with eigen-pairs
$\\{\lambda_{\ell},\phi_{\ell}(x)\\}_{\ell=1}^{\infty}$ such that
${\mathcal{C}}\phi_{\ell}(x)=\phi_{\ell}(x)\lambda_{\ell}$,
$\|\phi_{\ell}\|_{2}=1$ for all $\ell\in{\mathbb{N}}$ and
$\sum_{\ell=1}^{\infty}\lambda_{\ell}<\infty$, then we have the following
series representation for $u(x)$:

$u(x)=\sum_{\ell=1}^{\infty}u_{\ell}\phi_{\ell}(x),\quad
u_{\ell}:=\int_{D}u(x)\phi_{\ell}(x)\overset{ind}{\sim}\mathrm{q}\\!-\\!\mathrm{ED}(0,\lambda_{\ell})$
(14)

where $\mathrm{E}[u_{\ell}]=0$ and
$\mathrm{Cov}(u_{\ell},u_{\ell^{\prime}})=\lambda_{\ell}\delta_{\ell\ell^{\prime}}$
with Dirac function $\delta_{\ell\ell^{\prime}}=1$ if $\ell=\ell^{\prime}$ and
$0$ otherwise.

###### Proof.

Note we can approximate $\phi_{\ell}(x)\in L^{2}(D)$ with simple functions
$\tilde{\phi}_{\ell}(x)=\sum_{i=1}^{d}k_{i}\chi_{D_{i}}(x)$ where $D_{i}$’s
are measurable subsets of $D$ and $\chi_{D_{i}}(x)=1$ if $x\in D_{i}$ and $0$
otherwise. By the linear combination property of elliptic distributions [c.f.
Theorem 2.6.3 in 8],
$\tilde{u}_{\ell}=\int_{D}u(x)\tilde{\phi}_{\ell}(x)dx\sim\mathrm{q}\\!-\\!\mathrm{ED}(0,c)$
with $c=\alpha_{d}^{-1}\mathrm{E}[\tilde{u}_{\ell}^{2}]$ to be determined.
Note
$\alpha_{d}=\frac{2^{\frac{2}{q}}\Gamma(\frac{d}{2}+\frac{2}{q})}{d\Gamma(\frac{d}{2})}d^{1-\frac{2}{q}}$
comes from Proposition 3.2 and the scaling ${\bf
u}^{*}=d^{\frac{1}{2}-\frac{1}{q}}{\bf u}$ in Definition 3. We have
$\alpha_{d}=\frac{\Gamma(\frac{d}{2}+\frac{2}{q})}{\Gamma(\frac{d}{2})}\left(\frac{2}{d}\right)^{\frac{2}{q}}\to
1$ as $d\to\infty$. Taking the limit $d\to\infty$, we have
$u_{\ell}=\int_{D}u(x)\phi_{\ell}(x)dx\sim\mathrm{q}\\!-\\!\mathrm{ED}(0,c)$.
In general, by the similar argument we have

$\displaystyle\mathrm{Cov}(u_{\ell},u_{\ell^{\prime}})$
$\displaystyle=\mathrm{E}[u_{\ell}u_{\ell^{\prime}}]$
$\displaystyle=\int_{D}\int_{D}\mathrm{E}[u(x)u(x^{\prime})]\phi_{\ell}(x)\phi_{\ell^{\prime}}(x^{\prime})dxdx^{\prime}$
$\displaystyle=\int_{D}\int_{D}{\mathcal{C}}(x,x^{\prime})\phi_{\ell}(x)\phi_{\ell^{\prime}}(x^{\prime})dxdx^{\prime}$
$\displaystyle=\int_{D}\lambda_{\ell}\phi_{\ell}(x^{\prime})\phi_{\ell^{\prime}}(x^{\prime})dx^{\prime}=\lambda_{\ell}\delta_{\ell\ell^{\prime}}$

Thus it completes the proof. ∎

######  Remark 2.

If we factor $\sqrt{\lambda_{\ell}}$ into $u_{\ell}$, we have the following
expansion for Q-EP more comparable to (1) for Besov:

$u(x)=\sum_{\ell=1}^{\infty}\lambda_{\ell}^{-\frac{1}{2}}u_{\ell}\phi_{\ell}(x),\quad
u_{\ell}\overset{iid}{\sim}\mathrm{q}\\!-\\!\mathrm{ED}(0,1)\propto\pi_{q}(\cdot)$
(15)

Because of the definition (1) in terms of expanded series, there is no
explicit formula for prediction using Besov process. By contrast, prediction
with Q-EP can be done naturally through the following conditional distribution
[3, 8].

######  Proposition 3.3.

If ${\bf u}=({\bf u}_{1},{\bf
u}_{2})\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}(\bm{\mu},{\bf C})$ with
$\bm{\mu}=\begin{bmatrix}\bm{\mu}_{1}\\\ \bm{\mu}_{2}\end{bmatrix}$ and ${\bf
C}=\begin{bmatrix}{\bf C}_{11}&{\bf C}_{12}\\\ {\bf C}_{21}&{\bf
C}_{22}\end{bmatrix}$, then we have the following conditional distribution

$\displaystyle{\bf u}_{1}|{\bf u}_{2}$
$\displaystyle\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}(\bm{\mu}_{1\cdot 2},{\bf
C}_{11\cdot 2}),$ $\displaystyle\bm{\mu}_{1\cdot 2}$
$\displaystyle=\bm{\mu}_{1}+{\bf C}_{12}{\bf C}_{22}^{-1}({\bf
u}_{2}-\bm{\mu}_{2}),{\bf C}_{11\cdot 2}={\bf C}_{11}-{\bf C}_{12}{\bf
C}_{22}^{-1}{\bf C}_{21}$

###### Proof.

This directly follows from [Corollary 5 of Theorem 5 in 3] or [Corollary 3 of
Theorem 2.6.6 in 8] for $\mathrm{q}\\!-\\!\mathrm{ED}_{d}(\bm{\mu},{\bf
C})=\mathrm{EC}_{d}(\bm{\mu},{\bf C},g)$ with $g$ as chosen in (8). ∎

######  Remark 3.

Based on Proposition 3.2, we know that Q-EP has same predicative mean as GP
does. But their predicative covariance formulae differ by a constant.

##  4 NUMERICAL EXPERIMENTS

In this section, we first show that the elliptic slice sampler [ESS 18] can be
directly used to generate posterior samples of $u$ for models with the Q-EP
prior. Then we compare GP, Besov and Q-EP in modeling time series and
reconstructing images. These two numerical experiments demonstrate that our
proposed Q-EP enables faster convergence in obtaining a better maximum a
posterior (MAP) estimate. What is more, the posterior samples provide
uncertainty quantification by the posterior standard deviation, which is not
available in the frequentist methods.

###  4.1 Inference: Elliptic Slice Sampler

ESS [18] is an Markov chain Monte Carlo (MCMC) algorithm designed to sample
from posterior of models with GP priors. It first generates a random point
${\bf z}_{1}$ from the GP prior ${\mathcal{N}}({\bf 0},{\bf C})$ and with the
current state ${\bf z}_{0}$ it creates an ellipsis: ${\bf
z}_{0}\cos\theta+{\bf z}_{1}\sin\theta$ for $\theta\in[0,2\pi]$. Then it
continues with slice sampling [19] for the likelihood on the generated
ellipsis. The key to the proof of validity is that any point on the ellipsis
follows the same Gaussian distribution ${\mathcal{N}}({\bf 0},{\bf C})$. We
will show that the similar fact exists for the $q$-exponential distribution
$\mathrm{q}\\!-\\!\mathrm{ED}_{d}(0,{\bf C})$ thus ESS can be applied to
models with Q-EP priors.

######  Proposition 4.1.

If ${\bf u}_{0},{\bf u}_{1}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}({\bf 0},{\bf
C})$, then for ${\bf u}(\theta)={\bf u}_{0}\cos\theta+{\bf u}_{1}\sin\theta$
we have

${\bf u}(\theta)\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}({\bf 0},{\bf
C}),\quad\forall\theta\in[0,2\pi]$ (16)

###### Proof.

Based on the consistency, we have the joint distribution

$\begin{bmatrix}{\bf u}_{0}\\\ {\bf
u}_{1}\end{bmatrix}\sim\mathrm{q}\\!-\\!\mathrm{ED}_{2d}\left({\bf
0},\begin{bmatrix}{\bf C}&{\bf 0}\\\ {\bf 0}&{\bf C}\end{bmatrix}\right)$

By the linear combination property of the elliptic distributions [12, 8], we
have

$\displaystyle{\bf u}(\theta)$ $\displaystyle=\begin{bmatrix}{\bf
I}_{d}\cos\theta&{\bf I}_{d}\sin\theta\end{bmatrix}\begin{bmatrix}{\bf
u}_{0}\\\ {\bf u}_{1}\end{bmatrix}$
$\displaystyle\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}\left({\bf
0},\begin{bmatrix}{\bf I}_{d}\cos\theta&{\bf
I}_{d}\sin\theta\end{bmatrix}\begin{bmatrix}{\bf C}&{\bf 0}\\\ {\bf 0}&{\bf
C}\end{bmatrix}\begin{bmatrix}{\bf I}_{d}\cos\theta\\\ {\bf
I}_{d}\sin\theta\end{bmatrix}\right)$
$\displaystyle\sim\mathrm{q}\\!-\\!\mathrm{ED}_{d}({\bf 0},{\bf C})$

This completes the proof. ∎

###  4.2 Time Series Modeling

(a) Time series with step jumps

(b) Time series with sharp turnings.

Figure 2: MAP estimates by GP (left), Besov (middle) and Q-EP (right) models.
Blue solid lines are true trajectories. Orange dots are actual realizations
(data points). Red dashed lines are MAP estimates.

We first consider two time series, one with step jumps and the other with
sharp turnings, whose true trajectories are

$u_{s}(t)=\begin{cases}1,&t\in[0,1]\\\ 0.5,&t\in(1,1.5]\\\ 2,&t\in(1.5,2]\\\
0,&otherwise\end{cases}$ (17a) $u_{a}(t)=\begin{cases}1.5t,&t\in[0,1]\\\
3.5-2t,&t\in(1,1.5]\\\ 3t-4,&t\in(1.5,2]\\\ 0,&otherwise\end{cases}$ (17b)

(a) Time series with step jumps.

(b) Time series with sharp turnings.

Figure 3: Negative posterior densities (left) and errors (right) as functions
of iterations.

We generate these two time series $\\{y_{i}\\}$ by adding Gaussian noises to
their true trajectories evaluated at $N=200$ evenly spaced points
$\\{t_{i}\\}$ in $[0,2]$, that is,

$y_{i}^{*}=u_{*}(t_{i})+\varepsilon_{i},\quad\varepsilon_{i}\overset{ind}{\sim}N(0,\sigma^{2}_{*}(t_{i})),\quad
i=1,\cdots,N,\;*=s,a.$

Let $\sigma_{s}/\|f_{s}\|=0.015\;for\;t_{i}\in[0,2]$ and
$\sigma_{a}/\|f_{a}\|=0.01\;if\;t_{i}\in[0,1];\;0.07\;if\;t_{i}\in(1,2]$. See
Figures 2(a) and 2(b) for their true trajectories (blue lines) and
realizations (orange points) respectively.

(a) Time series with step jumps

(b) Time series with sharp turnings.

Figure 4: Predictions by GP (left) and Q-EP (right) models. Blue solid lines
are true trajectories. Black ticks indicate the training data points. Red
dashed lines are model estimates. Red dot-dashed lines are model predictions
with shaded region being credible bands.

We use the above likelihood and test three priors: GP, Besov and Q-EP. For
Besov, we choose the Fourier basis

$\phi_{0}(t)=\sqrt{2},\quad\phi_{\ell}(t)=\cos(\pi\ell
t),\quad\ell\in{\mathbb{N}}.$

For both GP and Q-EP, we adopt the following Matérn kernel with
$\nu=\frac{1}{2}$, $\sigma^{2}=1$, $l=0.5$ and $s=1$:

$\displaystyle C(t,t^{\prime})=$
$\displaystyle\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left[\sqrt{2\nu}(\|t-t^{\prime}\|/l)^{s}\right]^{\nu}\cdot$
$\displaystyle K_{\nu}(\sqrt{2\nu}(\|t-t^{\prime}\|/l)^{s})$

In both Besov and Q-EP, we set $q=1$. Figure 2(a) compares the MAP estimates
by GP, Besov and Q-EP for the time series with step jumps and Figure 2(a)
compares them for the time series with sharp turnings. We can see in both
cases, Q-EP yields the best estimates closest to the true trajectories. We
also investigate the negative posterior densities and relative errors,
$\|\hat{u}_{*}-u_{*}\|/\|u_{*}\|$, as functions of iterations in Figure 3. The
absolute values of negative posterior densities are not comparable since they
are different models. But their plots inform that both GP and Q-EP model
converge faster than Besov. The error reducing plots on the right panels in
both Figures 3(a) and 3(b) indicate that Q-EP prior model can achieve smaller
errors relative to the truth $u_{*}$ for $*=s,a$, verifying the better MAP
estimates in the previous figure 2.

Figure 5: Satellite image: true image (left) and the contaminated observation
(right).

Then we consider the prediction problem by holding out about $30\%$ of data
for testing. These testing data are selected by including the last $1/8$
portion and sub-sampling every other from the last but $3/8$ portion of the
whole data. The remaining $138$ data points, as indicated by short “ticks" in
Figure 4, are used to train the models with GP and Q-EP priors. This poses the
prediction model both interpolation (among observations) and extrapolation (at
no-observation region) tasks. We plot the prediction results by GP (left) and
Q-EP (right) in Figure 4. For both models, extrapolation (in the last $1/8$
portion of the data) comes with larger uncertainty (wider credible band) than
interpolation (in the last but $3/8$ portion of the data). Note, such
prediction is not immediately available for models based on the Besov prior.

###  4.3 Image Reconstruction

Figure 6: Satellite image: true image and MAP estimates by GP, Besov and Q-EP
models.

Next we consider reconstructing a ($128\times 128$ pixels) satellite image
shown on the left panel of Figure 5 from a contaminated observation on the
right panel. The image itself can be viewed as a function $u(x)$ on the square
unit $D=[0,1]^{2}$ taking values as the pixels. When evaluating $u(x)$ on the
discretized domain, $u(x)$ becomes a matrix of size $128\times 128$, which can
further be vectorized to ${\bf u}\in{\mathbb{R}}^{d}$ with $d=128^{2}$. The
true image, denoted as $u^{\dagger}$, is blurred by applying a motion blur
point spread function [PSF 2] and adding $5\%$ Gaussian noise. The actual
observation, $y(x)$, can be written as in the following linear model:

$y(x)=Au(x)+\varepsilon,\quad\varepsilon\sim{\mathcal{N}}(0,\sigma^{2}_{\varepsilon})$

where $A\in{\mathbb{R}}^{J\times d}$ is the blur motion PSF with $J=d$ and
$\sigma_{\varepsilon}/\|Au\|=5\%$.

In general, this linear inverse problem could be severely under-determined
($J\ll d$) thus it may become very challenging to reconstruct the true
satellite image $u^{\dagger}(x)$. Any reconstruction without proper
regularization, e.g. least square estimation, or ridge regression, will be far
away from being satisfactory.

Figure 7: Satellite image: negative posterior densities (left) and errors
(right) as functions of iterations.

We adopt the Bayesian approach to this inverse problem [6] and endow GP, Besov
and Q-EP as priors for the image function $u(x)$ respectively. Note, the
computation involving a full sized ($d\times d$) kernel matrix ${\bf C}$ for
GP and Q-EP is prohibitive. Therefore, we consider its Mercer’s expansion (13)
for a truncation with the first $L=2000$ items. We compare the reconstructions
by MAP estimate in Figure 6. The output by GP is too blurry and very close to
the observed image, which means that GP does not “de-noise" much from the
observation. The result by Besov is much better than GP due to the $L_{1}$
regularization. However, it is still a noisy reconstruction by Besov. We can
visually tell that the Q-EP prior model produces the reconstruction of the
highest quality. We also compare their negative posterior densities and
relative errors, $\|\hat{u}-u^{\dagger}\|/\|u^{\dagger}\|$, in Figure 7. Again
the Q-EP prior model yields the smallest error among the three models.

Table 1: Posterior estimates by samples $\\{u\\}$ of satellite image by GP, Besov and Q-EP. $\textrm{REM}:=\|\overline{u}-u^{\dagger}\|/\|u^{\dagger}\|$. Results are repeated for 10 times. | GP | Besov | Q-EP
---|---|---|---
REM | 0.3213 | 0.9071 | 0.3198
Std(REM) | 4.17e-6 | 1.54E-3 | 3.32E-3

Lastly, we apply ESS to these three models to obtain posterior samples of $u$.
Here we adopt the same white noise representation of Besov by [4] thus we can
also apply ESS to Besov prior models. We compare ESS with the whitened
preconditioned Crank-Nicolson [wpCN, proposed in 4] and find ESS outperforms
wpCN in obtaining better posterior reconstruction. Therefore we only include
ESS in the following comparisons. We run 10000 samples after burning 5000
samples and use them to reconstruct $u$ (posterior mean or median) and
quantify uncertainty (posterior standard deviation). Table 1 summarizes the
relative error of mean,
$\textrm{REM}:=\|\overline{u}-u^{\dagger}\|/\|u^{\dagger}\|$ (with
$\overline{u}$ being the mean of posterior samples), and the standard
deviations of REM’s for repeating the experiments 10 times by three prior
models respectively. Q-EP attains the lowest relative error $0.3198$.

Figure 8: Satellite image: uncertainty field (posterior standard deviation)
given by GP, Besov and Q-EP models. They not plotted in the same scale: GP has
much smaller uncertainty values (about
$1{}^{\text{o}}\mkern-5.0mu/\mkern-3.0mu_{\text{oo}}$ of that with Q-GP).

In Figure 8, we compare the uncertainty field (posterior standard deviation)
by these three prior models. It appears that GP has more recognizable
uncertainty filed than the other two. However, they are not plotted in the
same scale: GP has much smaller posterior standard deviation (about
$1{}^{\text{o}}\mkern-5.0mu/\mkern-3.0mu_{\text{oo}}$ of that with Q-GP)
compared with the other two. Therefore, this raises a red flag that GP could
be over-confident about a less accurate estimate.

##  5 CONCLUSION

In this paper we propose the $q$-exponential process (Q-EP) as a prior on
$L^{q}$ functions with a flexible parameter $q>0$ to control the degree of
regularization. Usually, $q=1$ is adopted to capture abrupt changes or sharp
contrast in data such as edges in the image as the Besov prior has recently
gained popularity for. Motivated by the same $q$-exponential distribution in
the definition of Besov process, Q-EP can be viewed as a probabilistic
formulation of Besov process with direct control on the correlation length and
an explicit formula to make prediction. Compared with GP, Q-EP can impose
sharper regularization through $q$. Compared with Besov, Q-EP enjoys the
explicit formula with more control on the correlation structure as GP. The
numerical experiments in time series modeling and image reconstruction
demonstrate our proposed Q-EP is superior in obtaining better reconstruction
faster.

In the current work, we manually fix hyper-parameters in the kernel. The
reported results are not sensitive to some of these hyper-parameters such as
the magnitude $\sigma^{2}$ and the correlation length $l$ but may change
drastically to others like the regularity parameter $\nu$ and the smoothness
parameter $s$. In future work, we will incorporate hyper-priors for some of
those parameters and adopt a hierarchical scheme. We plan to study the
properties such as regularity of function draws of Q-EP and the posterior
contraction. Future work will also consider operator based kernels such as
graph Laplacian [6, 7, 15].

### Acknowledgement

SL is supported by NSF grant DMS-2134256. We thank Mirjeta Pasha for sharing
the data of satellite image.

## References

  * [1] Maria Bånkestad, Jens Sjölund, Jalil Taghia, and Thomas Schön.  The elliptical processes: a family of fat-tailed stochastic processes.  03 2020.
  * [2] Alessandro Buccini, Mirjeta Pasha, and Lothar Reichel.  Linearized krylov subspace bregman iteration with nonnegativity constraint.  Numerical Algorithms, 87(3):1177–1200, sep 2020.
  * [3] Stamatis Cambanis, Steel Huang, and Gordon Simons.  On the theory of elliptically contoured distributions.  Journal of Multivariate Analysis, 11(3):368–385, 1981.
  * [4] Victor Chen, Matthew M. Dunlop, Omiros Papaspiliopoulos, and Andrew M. Stuart.  Dimension-robust mcmc in bayesian inverse problems.  03 2018.
  * [5] Masoumeh Dashti, Stephen Harris, and Andrew Stuart.  Besov priors for bayesian inverse problems.  Inverse Problems and Imaging, 6(2):183–200, may 2012.
  * [6] Masoumeh Dashti and Andrew M. Stuart.  The Bayesian Approach to Inverse Problems, pages 311–428.  Springer International Publishing, Cham, 2017.
  * [7] Matthew M. Dunlop, Dejan Slepčev, Andrew M. Stuart, and Matthew Thorpe.  Large data and zero noise limits of graph-based semi-supervised learning algorithms.  Applied and Computational Harmonic Analysis, 49(2):655–697, 2020\.
  * [8] K. Fang and Y.T. Zhang.  Generalized Multivariate Analysis.  Science Press, 1990.
  * [9] E. Gómez, M.A. Gomez-Viilegas, and J.M. Marín.  A multivariate generalization of the power exponential family of distributions.  Communications in Statistics - Theory and Methods, 27(3):589–600, jan 1998.
  * [10] A. P. Dawid J. M. Bernardo, J. O. Berger and A. F. M. Smith.  Regression and classification using gaussian process priors.  Bayesian Statistics, 6:475–501, 1998.
  * [11] Junxiong Jia, Jigen Peng, and Jinghuai Gao.  Bayesian approach to inverse problems for functions with a variable-index besov prior.  Inverse Problems, 32(8):085006, 2016.
  * [12] Mark E. Johnson.  Multivariate Statistical Simulation, chapter 6 Elliptically Contoured Distributions, pages 106–124.  Probability and Statistics. John Wiley & Sons, Ltd, 1987.
  * [13] Y. Kano.  Consistency property of elliptic probability density functions.  Journal of Multivariate Analysis, 51(1):139–147, 1994.
  * [14] Tomasz J. Kozubowski, Krzysztof Podgórski, and Igor Rychlik.  Multivariate generalized laplace distribution and related random fields.  Journal of Multivariate Analysis, 113:59–72, 2013.  Special Issue on Multivariate Distribution Theory in Memory of Samuel Kotz.
  * [15] Shiwei Lan.  Learning temporal evolution of spatial dependence with generalized spatiotemporal gaussian process models.  Journal of Machine Learning Research, 23(259):1–53, 2022.
  * [16] Matti Lassas, Eero Saksman, and Samuli Siltanen.  Discretization-invariant bayesian inversion and besov space priors.  Inverse Problems and Imaging, 3(1):87–122, 2009.
  * [17] Felix Lucka.  Fast markov chain monte carlo sampling for sparse bayesian inference in high-dimensional inverse problems using l1-type priors.  Inverse Problems, 28(12):125012, nov 2012.
  * [18] Iain Murray, Ryan Prescott Adams, and David J. C. MacKay.  Elliptical slice sampling.  9:541–548, 2010.
  * [19] Radford M. Neal.  Slice sampling.  The Annals of Statistics, 31(3), jun 2003.
  * [20] Bernt Øksendal.  Stochastic Differential Equations.  Springer Berlin Heidelberg, 2003.
  * [21] Krzysztof Podgórski and Jörg Wegener.  Estimation for stochastic models driven by laplace motion.  Communications in Statistics - Theory and Methods, 40(18):3281–3302, sep 2011.
  * [22] Carl Edward Rasmussen and Christopher K. I. Williams.  Gaussian Processes for Machine Learning.  The MIT Press, 2005.
  * [23] I. J. Schoenberg.  Metric spaces and completely monotone functions.  Annals of Mathematics, 39:811–841, 1938.
  * [24] Simopekka Vänskä, Matti Lassas, Samuli Siltanen, and Rolf Insitute.  Statistical x-ray tomography using empirical besov priors.  International Journal of Tomography and Statistics, 11, 06 2009\.


# An Orthogonal Polynomial Kernel-Based Machine Learning Model for
Differential-Algebraic Equations

Tayebeh Taheri Email: {ttaherii1401<EMAIL_ADDRESS>Department
of Computer and Data Sciences, Faculty of Mathematical Sciences, Shahid
Beheshti University, G.C. Tehran, Iran Alireza Afzal Aghaei††footnotemark:
Department of Computer and Data Sciences, Faculty of Mathematical Sciences,
Shahid Beheshti University, G.C. Tehran, Iran Kourosh Parand Email:
<EMAIL_ADDRESS>Corresponding author Department of Computer and Data
Sciences, Faculty of Mathematical Sciences, Shahid Beheshti University, G.C.
Tehran, Iran Department of Cognitive Modeling, Institute for Cognitive and
Brain Sciences, Shahid Beheshti University, G.C. Tehran, Iran  Department of
Statistics and Actuarial Science, University of Waterloo, Waterloo, Canada

###### Abstract

The recent introduction of the Least-Squares Support Vector Regression (LS-
SVR) algorithm for solving differential and integral equations has sparked
interest. In this study, we expand the application of this algorithm to
address systems of differential-algebraic equations (DAEs). Our work presents
a novel approach to solving general DAEs in an operator format by establishing
connections between the LS-SVR machine learning model, weighted residual
methods, and Legendre orthogonal polynomials. To assess the effectiveness of
our proposed method, we conduct simulations involving various DAE scenarios,
such as nonlinear systems, fractional-order derivatives, integro-differential,
and partial DAEs. Finally, we carry out comparisons between our proposed
method and currently established state-of-the-art approaches, demonstrating
its reliability and effectiveness.

Keywords— Least-squares support vector regression, Differential-algebraic
equations, Fractional derivative, Machine learning methods, Weighted residual
methods

##  1 Introduction

A system of differential-algebraic equations (DAEs) is a combination of
differential equations and algebraic equations, in which the differential
equations are related to the dynamical evolution of the system, and the
algebraic equations are responsible for constraining the solutions that
satisfy the differential and algebraic equations. DAEs serve as essential
models for a wide array of physical phenomena. They find applications across
various domains such as mechanical systems, electrical circuit simulations,
chemical process modeling, dynamic system control, biological simulations, and
control systems. Consequently, solving these intricate differential equations
has remained a significant challenge for researchers. To address this, a range
of techniques including numerical, analytical, and semi-analytical methods
have been employed to tackle the complexities inherent in solving DAEs.
Examples of these techniques include Adomian Decomposition Method (ADM) [7],
Pade series [6], Runge-Kutta [4], variational iteration [18], homotopy
perturbation [19], power series method [3], homotopy analysis method [23],
spectral methods [2], and so on. Table 1 presents a comparison of properties
among spectral, deep learning, and LS-SVR methods.

In recent years, machine learning schemes have been developed to solve various
differential equations, such as ordinary differential equations (ODEs),
partial differential equations (PDEs), stochastic differential equations, and
more. These innovative approaches try to extract the problem’s underlying
structure by learning from the given differential equations. The support
vector machine technique is one of the most powerful machine learning methods
that, according to a set of labeled training examples, tries to increase the
generalization of the learned solution. Many scientists and researchers have
used this method to solve differential equations. For instance, Hajimohammadi
et al. [12] proposed a numerical learning method to solve the general
Falkner–Skan model. A least-squares support vector regression scheme has been
developed for solving differential and integral equations by Aghaei and Parand
[1]. Lu et al. [14] solved a higher nonlinear ordinary differential equation
via a method based on least-squares support vector machines. In the mentioned
articles, we see instances of the use of support vector regression to solve
different types of differential equations. However, until now, this method
[21] has not been used to solve a system of differential-algebraic equations.
According to the noted statements, in this research, we seek to develop this
methodology for solving systems of differential-algebraic equations. A
significant point of consideration within this paper is that these equations
have been solved with fractional derivative orders, the solutions of which are
considerably more intricate in comparison to equations involving integer
derivatives. Since a system of equations combines both differential and
algebraic components, the process of resolving DAE systems inherently entails
a higher level of complexity when compared to the task of solving ODE systems.
In this paper, the application of weighted residual methods (WRM) to
approximate solutions of differential equations is utilized. A cohesive
integration of the collocation approach and the least-squares support vector
regression method has been employed to determine the weight coefficients.
Notably, within this framework, the adoption of Legendre functions as kernel
functions in the context of the least squares support vector regression
methodology is a distinctive facet. The proposed methodology has been employed
to achieve numerical solutions for various categories of systems of DAEs,
including fractional DAEs, integro-differential algebraic equations of
fractional order, and partial differential-algebraic equations.

In brief, the distinguishing features of our work encompass the following
points:

  * •

Merging the collocation method with the LS-SVR approach.

  * •

Utilizing Legendre functions as kernels within the LS-SVR algorithm.

  * •

Utilizing LS-SVR for addressing linear and nonlinear DAEs involving fractional
derivatives.

  * •

Utilizing LS-SVR for addressing Volterra integro-differential algebraic
equations of fractional order.

  * •

Utilizing LS-SVR for addressing partial DAEs.

Property |  Spectral Methods |  Deep Learning |  LS-SVR
---|---|---|---
Convergence |  Exponential Convergence. |  Ensuring convergence to the best solution is not always guaranteed. |  Exponential Convergence.
Time Complexity |  Solving a linear system of equations. |  Usually computationally expensive, especially for large models. |  Solving a positive definite linear system of equations.
Parameters |  Few parameters. |  Parameters, including weights and biases, can be numerous |  Few parameters.
Scalability |  Calculations become challenging when dealing with differential equations in high dimensions. |  Scalable for high dimensions and complex domains. |  Calculations are superior to spectral methods but fall short of the performance achieved by deep learning.
Table 1: Comparison of Methods

The rest of this paper is organized as follows. In section 2, first, some
types of differential-algebraic equations are reviewed, and a brief
explanation of fractional derivative calculations is given. Section 3
introduces the least-squares support vector machine regression model,
following a presentation of orthogonal polynomials and their characteristics.
The efficacy of our approach is demonstrated through its implementation in
approximating various examples, as detailed in section 4. Finally, a summary
of the proposed method is given in section 5.

##  2 Preliminaries

This section explains some algebraic-differential equations solved in the
numerical examples section. Ultimately, we introduce Caputo’s fractional
derivative for differential equations with fractional derivative.

###  2.1 A system of differential-algebraic equations

A system of differential-algebraic equations is a mathematical model where
certain variables are described by differential equations, while others are
constrained by algebraic equations.

The general form of this case is:

$\displaystyle f(t,X(t),X^{\prime}(t),Y(t))=0,$ (1) $\displaystyle
g(t,X(t),Y(t))=0,$ (2)

where $t$ is independent variable, $X\in\mathbb{R}^{n}$ is a vector of
differential variables, $X^{\prime}$ is derivative of $X$,
$Y\in\mathbb{R}^{m}$ is the vector of algebraic variables, function $f$ is a
set of algebraic equations, and function $g$ describes the algebraic
constraints of the system. These equations are known based on the index of the
DAE. The index is the highest derivative that the differential variables
should be differentiated. The form of the DAE is also classified based on the
type of algebraic equations it includes, such as semi-explicit and fully
implicit DAEs. In Section 4, we address the resolution of an index-1 semi-
explicit DAE [13] within Example 1. The general form of these kinds of
differential equations is given by:

$\begin{cases}X^{\prime}(t)=f(t,X(t),Y(t))\\\ g(t,X(t),Y(t))=0.\\\
\end{cases}.$ (3)

###  2.2 Integro-differential algebraic equations

Integro-differential algebraic equations (IDAEs) combine differential
equations and algebraic equations with integrals. In this system, some
equations have integrals of the dependent variables. The general form of this
case can be expressed as [22]:

$f(t,X^{(v)}(t),Y(t))+\int_{a}^{t}h(t,s,X^{(\leq
v)},Y(s))ds=g(t,X(t),Y(t))=0,$ (4)

where $t\in[a,t_{f}]$ is the independent variable and with dependent variable
$X$ which their derivatives up to $v$-th order. The equation (4) is made of
one section of differential-algebraic equation, $f(t,X^{(v)}(t),Y(t))$, and
one section of integral algebraic equation, $h(t,s,X^{(v)},Y(s)$.

  * •

If $h=0$, the equation (4) will be a system of differential–algebraic
equations.

  * •

If $f=0$, the equation (4) will be a integral algebraic equation.

  * •

If $v=0$, the equation (4) will be a system of integral–algebraic equations.

In section 4, within example 2, a case study focuses on solving a Volterra
integro-differential algebraic equation (VI-DAE) with weakly singular kernels
of index-1.

###  2.3 Partial differential-algebraic equations

Partial differential-algebraic equations (PDAEs) are a class of mathematical
equations that involve one partial differential equation, one ordinary
differential equation, and one algebraic equation. Partial differential
equations have partial derivatives related to the spatial coordinates.
However, PDAEs have derivatives concerning both spatial and time coordinates.
This makes PDAEs much more complex than PDEs and requires a different
approach. This is a form of linear partial differential-algebraic equations
[11]:

$Av_{t}(t,x)+Bv_{xx}(t,x)+Cv(t,x)=g(t,x),$ (5)

where $A,B,C\in\mathbb{R}^{n\times xn}$ are constant matrices that can be
singular, $t\in(0,t_{e}),x\in(-l,l)\subset\mathbb{R}$, and
$v,g:[0,t_{e}]\times[-l,l]\rightarrow\mathbb{R}^{n}$.

  * •

If $A=0$, the equation (5) will be an ordinary differential equation.

  * •

If $B=0$, the equation (5) will be a differential-Algebraic equation.

In section 4, as illustrated in example 5, a partial differential-algebraic
equation has been approximated. In the next part DAEs with fractional order
will be explained.

###  2.4 Fractional differential-algebraic equations

Fractional differential-algebraic equations (FDAEs) involve fractional
derivatives, which can lead to precise and adaptable for specific physical
systems to anticipate and regulate their behavior. Section 4 covers diverse
categories of FDAEs, encompassing linear, nonlinear, and Volterra integro-
differential algebraic equations. In general, an FDAE can be written in the
following form:

$\displaystyle\prescript{}{}{\mathcal{D}}^{\alpha}x(t)=f(t,X(t),D^{v}X(t),Y(t))$
(6) $\displaystyle g(t,X(t),D^{v}X(t),Y(t))=0,$ (7)

where $\prescript{}{}{\mathcal{D}}^{\alpha}$ denotes the fractional derivative
of order $\alpha$ which may be given in terms of Caputo or Riemann-Liouville
fractional derivative. In the next part, a brief explanation of calculating
the fractional derivative is provided.

######  Definition 1.

Consider a continuous function $h$ over the interval $[a,b]$. The Riemann-
Liouville fractional derivative of any order $\alpha$ is defined as follows:

$\prescript{}{}{\mathcal{D}}^{\alpha}h(x)=\frac{1}{\Gamma(m-\beta)}\frac{d^{m}}{dx^{m}}\int_{a}^{x}(x-t)^{m-\alpha-1}h(t)dt,\quad
m-1\leq\alpha\leq m,$ (8)

here $\Gamma(.)$ represents the Gamma function, symbolized by the subsequent
integral:

$\Gamma(z)=\int_{0}^{\infty}e^{-\tau}\tau^{z-1}d\tau,$ (9)

that this integral is convergent in the right half of the complex plane
$(\text{Re}(z)\geq 0)$.

Caputo’s [5] formulation of non-integer order differentiation presents a
distinctive approach to the concept. A notable advantage of this definition is
its ability to express initial conditions of fractional differential equations
in a format consistent with those for integer-order differential equations.
The following definition offers an alternative approach to formulating the
fractional derivative.

######  Definition 2.

Caputo’s formulation for the fractional order $\alpha$ is presented as:

$\prescript{C}{}{\mathcal{D}}^{\alpha}h(x)=\frac{1}{\Gamma(v-\alpha)}\int_{a}^{x}\frac{h^{(v)}(t)dt}{(x-t)^{\alpha+1-v}},\quad
v-1<\alpha<v.$ (10)

In Caputo’s derivative formula, the calculation of this derivative requires
the calculation of an integral, which may be complicated and computationally
inefficient. Therefore, in the next theorem, a method to approximate this
derivative is introduced.

######  Theorem 1.

Let $0<\alpha<1$ and the interval $[x_{0},x_{m}]$ is discretized to $m$
points, $0=x_{0}<x_{1}<\dots<x_{m}$. By subtracting an integral part in
Caputo’s definition to $m$ sub-interval, Caputo’s fractional derivative is
defined as:

$\prescript{C}{}{\mathcal{D}}^{\alpha}h(x_{m})=\frac{1}{\Gamma(m-\alpha)}\sum_{k=0}^{m-1}\int_{x_{k}}^{x_{k+1}}\frac{h^{\prime}(t)}{({x_{m}-t})^{\alpha}}dt.$
(11)

And since the first derivative of the function $h(x)$ is estimated as:

$h^{\prime}(t)\approx\frac{h_{k+1}-h_{k}}{f_{k+1}-f_{k}},\quad
x_{k}<t<x_{k+1},$ (12)

that $h_{k}=h(f_{k})$, then the formula (11) can be expressed as:

$\prescript{C}{}{\mathcal{D}}^{\alpha}h(x_{m})\approx\sum_{k=0}^{m}(g_{k-1}-g_{k})h_{k},$
(13)

with $g_{k}$ defined as:

$g_{k}=\begin{cases}\frac{(x_{m}-x_{k})^{1-\alpha}-(x_{m}-x_{k+1})^{1-\alpha}}{\Gamma(2-\alpha)(x_{k+1}-x_{k}}),&1\leq
k<m\\\ 0,&otherwise\end{cases}$.

This process is known as $L1$ discretization, in which $\alpha$ is the order
of the fractional derivative and assumes that the points $\Delta
x=x_{k+1}-x_{k}$ are equidistant.

###### Proof.

[15] ∎

Fractional calculus has led to the development of different mathematical
models, which have better consistency with experimental data in both physical
and engineering problems, like fractional differential equations.

##  3 Methodology

This section proposes a model for solving differential-algebraic equations
using the least-squares support vector regression method in conjunction with
the collocation approach and the Legendre polynomials.

###  3.1 Legendre Orthogonal Kernel

As a weighted residual approach, Spectral methods provide the advantage of
choosing basis functions. In spectral methods, basis functions are expected to
be infinitely differentiable and possess a completeness property. Commonly,
orthogonal functions are employed as the basis, simplifying calculations and
enabling the derivation of structured systems. Within this context, Jacobi
polynomials, a family of orthogonal polynomials defined over the interval
$[-1,1]$ with respect to the weight function $w(x)$, include the renowned
Legendre polynomial sequence. The Legendre polynomials serve as global
functions, constituting a comprehensive basis for square-integrable functions
on the interval $[-1,1]$. Their utilization leads to spectral convergence,
enhancing solution accuracy as the degree of Legendre polynomials rises.
Furthermore, explicit Legendre formulas enable efficient computation in
numerical methods, making them an ideal choice for solving differential
equations.

######  Definition 3.

Two polynomials from an infinite sequence of polynomials
$\left\\{P_{n}(x)\right\\}^{\infty}_{n=0}$ are considered orthogonal if their
inner product, determined for an integrable weight function $w(x)\geq 0$ over
the interval $[a,b]$, evaluates to zero. Furthermore, Legendre polynomials
associated with the weight function $w(x)=1$ within the interval $[-1,1]$ are
defined by the formula [9]:

$\langle{P_{n}(x),P_{m}(x)}\rangle=\frac{2}{2n+1}\delta_{nm},$ (14)

where $\delta$ signifies the Kronecker delta function. The following formula
presents the definition of the Legendre polynomial with degree $n$.

######  Definition 4.

These polynomials can be obtained through the expression [17]:

$P_{n}(x)=\sum_{v=0}^{\left[\frac{n}{2}\right]}(-1)^{v}\frac{(2n-2v)!}{2^{n}(n-v)!(n-2v)!v!}t^{n-2v},$
(15)

where $v$ and $n$ are non-negative integers. For enhanced clarity, we
calculate Legendre polynomials of degrees zero to three using (15):

$\displaystyle P_{0}(x)=1$ $\displaystyle P_{1}(x)=x$ $\displaystyle
P_{2}(x)=\frac{1}{2}(3x^{2}-1)$ $\displaystyle
P_{3}(x)=\frac{1}{2}(5x^{3}-3x).$

Legendre polynomials can be transferred to arbitrary intervals with
appropriate mappings that preserve the orthogonality property. In order to use
shifted Legendre polynomials in the interval $[a,b]$, it is necessary to use
an affine transformation. For this purpose, the linear mapping

$\phi(t)=\frac{2t-a-b}{b-a},$ (16)

that maps the interval $[-1,1]$ to the interval $[a,b]$ is applied.

###  3.2 LS-SVR for solving differential-algebraic equations

In this section, we are trying to solve a system of differential-algebraic
equations by using the least-squares support vector regression model and
combining it with weighted residual methods. First, we begin by recalling the
LS-SVR model in the context of regression problems.

######  Theorem 2.

Consider the training data $\left\\{(x_{k},y_{k})\right\\}^{N}_{k=1}$, where
$x_{k}\in\mathbb{R}^{d}$ and $y\in\mathbb{R}$ are assumed. Within the
framework of the least-squares support vector machine, the unknown function is
formulated as follows:

$y(x)=\mathbf{w}^{T}\phi(x)+b,$ (17)

where $\mathbf{w}$ represents a weight vector, and $b$ corresponds to a bias
term. The primal form of the optimization problem is expressed as follows:

$\displaystyle\min\limits_{\mathbf{w},e,b}\quad$
$\displaystyle\frac{1}{2}\mathbf{w}^{T}\mathbf{w}+\gamma\frac{1}{2}\sum_{k=1}^{N}e_{k}^{2}$
s.t. $\displaystyle y_{k}=\mathbf{w}^{T}\phi(x)+b+e_{k},\quad k=1,\ldots,N,$

where, $e_{k}$ denotes the error variable, and $\gamma$ signifies the
regularization parameter. Through the utilization of Lagrange multipliers and
optimality conditions, the resultant system of linear equations for solving
the regression problem is presented below:

$\displaystyle\left[\begin{array}[]{c|c}0&1_{N}^{T}\\\ \hline\cr
1_{N}&\Omega+I/\gamma\end{array}\right]\left[\begin{array}[]{c}b\\\
\boldsymbol{\alpha}\end{array}\right]=\left[\begin{array}[]{c}0\\\
y\end{array}\right],$

here, $\Omega_{kl}=K(x_{k},x_{l})$ for all $k,l=1,\dots,N$, and
$\boldsymbol{\alpha}=\left[\alpha_{1},\alpha_{2},\dots,\alpha_{N}\right]^{T}$
represents the Lagrange multipliers. Also, we have:

$y(x)=\sum_{k=1}^{N}\alpha_{k}K(x_{k},x_{l})+b,$ (18)

which the function $K(x_{k},x_{l})$ is the kernel function.

###### Proof.

[20] ∎

To solve the system of differential-algebraic equations, a set of $k$
equations is characterized by corresponding unknown functions $f_{i}$:

$\mathcal{N}_{i}(u_{1},u_{2},\ldots,u_{i})=f_{i},\quad i=1,2,\ldots,k.$ (19)

The approximated solutions of $u_{i}$ can be expressed as a formulation
derived from the research of Suykens et al. [20]:

$u_{i}\approx\tilde{u}_{i}(x)=w_{i}^{T}\varphi(x)+b_{i}=\sum_{j=0}^{k}w_{j}\varphi_{j}(x)+b_{j},\quad
j=1,2,\ldots,k,$

where $\varphi_{i}(x)$ is replaced with Legendre functions of degree $i$, and
$w_{i}$ and $b_{i}$ are the unknown coefficients. In the next step, the
unknown coefficients are set side by side in a vector form. Thus, we have:

${\mathbf{w}}=[w_{1,1},w_{1,2},\ldots,w_{1,d},w_{2,1},w_{2,2},\ldots,w_{2,d},\ldots,w_{k,1},w_{k,2},\ldots,w_{k,d}],$

as observed within $d$-dimensional vectors. In order to find the unknown
coefficients $w_{i}$ and $b_{i}$ the optimization problem by using the inner
product, which is defined in (14), for solving (19) can be constructed as:

$\displaystyle\min\limits_{\mathbf{w},\mathbf{e}}$
$\displaystyle\quad\frac{1}{2}\mathbf{w}^{T}\mathbf{w}+\frac{\gamma}{2}\mathbf{e}^{T}\mathbf{e}$
(20) s.t.
$\displaystyle\quad\langle\mathcal{N}_{i}(\tilde{u}_{1},\tilde{u}_{2},\ldots,\tilde{u}_{k})-f_{i},\psi_{j}\rangle=e_{i,j},$
$\displaystyle j=1,\ldots,n,$

where $i=1,2,\ldots,k$ shows that $kn$ is the number of training points,
$e_{i,j}$ is the residual value, and $\psi_{j}$ is a set of basis functions in
the test space. The matrix $e_{i,j}$ should be vectorized the same as unknown
coefficients $\mathbf{w}$. For any linear operator $\mathcal{N}$, denoted by
$\mathcal{L}$, the dual form of the optimization problem (20) can be derived.
Here we obtain the dual form for a system of three equations and three unknown
functions. This process can be generalized for an arbitrary number of
equations. Suppose the following system of equations is given:

$\begin{cases}\mathcal{L}_{1}(u_{1},u_{2},u_{3})=f_{1}\\\
\mathcal{L}_{2}(u_{1},u_{2},u_{3})=f_{2}\\\
\mathcal{L}_{3}(u_{1},u_{2},u_{3})=f_{3}\end{cases}.$ (21)

By approximating the solutions using,

$\displaystyle\tilde{u}_{1}(x)$ $\displaystyle=w_{1}^{T}\varphi(x)+b_{1}.$
(22) $\displaystyle\tilde{u}_{2}(x)$
$\displaystyle=w_{2}^{T}\varphi(x)+b_{2}.$ $\displaystyle\tilde{u}_{3}(x)$
$\displaystyle=w_{3}^{T}\varphi(x)+b_{3}.$

Note that, for partial differential-algebraic equations, $\tilde{u}_{i}$ is
defined as $\tilde{u}_{i}(x,t)=w_{i,j}{\varphi}^{T}_{j}(x)\varphi_{i}(t)$,
which applies to the general form we present. The optimization problem (20)
takes the form

$\displaystyle\min\limits_{\mathbf{w},\mathbf{e}}$
$\displaystyle\quad\frac{1}{2}\mathbf{w}^{T}\mathbf{w}+\frac{\gamma}{2}\mathbf{e}^{T}\mathbf{e}$
s.t.
$\displaystyle\quad\langle\mathcal{L}_{1}(\tilde{u}_{1},\tilde{u}_{2},\tilde{u}_{3})-f_{1},\psi_{j}\rangle=e_{1,j},$
$\displaystyle j$ $\displaystyle=1,\ldots,n,$ s.t.
$\displaystyle\quad\langle\mathcal{L}_{2}(\tilde{u}_{1},\tilde{u}_{2},\tilde{u}_{3})-f_{2},\psi_{j}\rangle=e_{2,j},$
$\displaystyle j$ $\displaystyle=n+1,\ldots,2n,$ s.t.
$\displaystyle\quad\langle\mathcal{L}_{3}(\tilde{u}_{1},\tilde{u}_{2},\tilde{u}_{3})-f_{3},\psi_{j}\rangle=e_{3,j},$
$\displaystyle j$ $\displaystyle=2n+1,\ldots,3n,$

where

$\displaystyle\mathbf{w}$
$\displaystyle=[w_{1},w_{2},w_{3}]=[w_{1,1},w_{1,2},\ldots,w_{1,n},w_{2,1},w_{2,2},\ldots,w_{2,n},w_{3,1},w_{3,2},\ldots,w_{3,n}]$
(23) $\displaystyle\mathbf{e}$
$\displaystyle=[e_{1},e_{2},e_{3}]=[e_{1,1},e_{1,2},\ldots,e_{1,n},e_{2,1},e_{2,2},\ldots,e_{2,n},e_{3,1},e_{3,2},\ldots,e_{3,n}].$

In order to solve this system, the Lagrangian function should be constructed
as follows:

$\displaystyle\mathfrak{L}(\mathbf{w},\mathbf{e},\boldsymbol{\alpha})$
$\displaystyle=\frac{1}{2}\mathbf{w}^{T}\mathbf{w}+\frac{\gamma}{2}\mathbf{e}^{T}\mathbf{e}-\sum_{j=1}^{n}\alpha_{j}\langle\mathcal{L}_{1}(\tilde{u}_{1},\tilde{u}_{2},\tilde{u}_{3})-f_{1},\psi_{j}\rangle-
e_{j}-\sum_{j=1}^{n}\alpha_{n+j}\langle\mathcal{L}_{2}(\tilde{u}_{1},\tilde{u}_{2},\tilde{u}_{3})-f_{2},\psi_{j}\rangle-
e_{n+j}$
$\displaystyle-\sum_{j=1}^{n}\alpha_{2n+j}\langle\mathcal{L}_{3}(\tilde{u}_{1},\tilde{u}_{2},\tilde{u}_{3})-f_{3},\psi_{j}\rangle-
e_{2n+j}.$

The conditions for optimality of the Lagrangian function are partial
derivatives $\frac{\partial\mathfrak{L}}{\partial
w_{k}},\frac{\partial\mathfrak{L}}{\partial
e_{k}},\frac{\partial\mathfrak{L}}{\partial b}$, and
$\frac{\partial\mathfrak{L}}{\partial\alpha_{k}}$ which should be equaled
zero. Then we compute these equations to find the coefficients:

$\displaystyle\frac{\partial\mathfrak{L}}{\partial w_{k}}=0\rightarrow
w_{k}=\begin{cases}\displaystyle\sum_{j=1}^{n}\alpha_{j}\langle\mathcal{L}_{1}(\varphi_{k},0,0)-f_{1},\psi_{j}\rangle-
e_{j}+\sum_{j=1}^{n}\alpha_{n+j}\langle\mathcal{L}_{2}(\varphi_{k},0,0)-f_{2},\psi_{j}\rangle-
e_{n+j}+\\\
\displaystyle\sum_{j=1}^{n}\alpha_{2n+j}\langle\mathcal{L}_{3}(\varphi_{k},0,0)-f_{3},\psi_{j}\rangle-
e_{2n+j},&k=1,2,\ldots,d.\\\
\displaystyle\sum_{j=1}^{n}\alpha_{j}\langle\mathcal{L}_{1}(0,\varphi_{k},0)-f_{1},\psi_{j}\rangle-
e_{j}+\sum_{j=1}^{n}\alpha_{n+j}\langle\mathcal{L}_{2}(0,\varphi_{k},0)-f_{2},\psi_{j}\rangle-
e_{n+j}+\\\
\displaystyle\sum_{j=1}^{n}\alpha_{2n+j}\langle\mathcal{L}_{3}(0,\varphi_{k},0)-f_{3},\psi_{j}\rangle-
e_{2n+j},&k=1,2,\ldots,d.\\\
\displaystyle\sum_{j=1}^{n}\alpha_{j}\langle\mathcal{L}_{1}(0,0,\varphi_{k})-f_{1},\psi_{j}\rangle-
e_{j}+\sum_{j=1}^{n}\alpha_{n+j}\langle\mathcal{L}_{2}(0,0,\varphi_{k})-f_{2},\psi_{j}\rangle-
e_{n+j}+\\\
\displaystyle\sum_{j=1}^{n}\alpha_{2n+j}\langle\mathcal{L}_{3}(0,0,\varphi_{k})-f_{3},\psi_{j}\rangle-
e_{2n+j},&k=1,2,\ldots,d.\\\ \end{cases}$ (24)
$\displaystyle\frac{\partial\mathfrak{L}}{\partial e_{k}}$
$\displaystyle=0\rightarrow\gamma e_{k}+\alpha_{k}=0,\qquad k=1\ldots 3n.$
(25) $\displaystyle\frac{\partial\mathfrak{L}}{\partial b}$
$\displaystyle=0\rightarrow\begin{cases}\displaystyle\frac{\partial\mathfrak{L}}{\partial
b_{1}}&=0\rightarrow\displaystyle\sum_{i=1}^{n}\alpha_{i}\langle\mathcal{L}_{1}(1,0,0),\psi_{i}\rangle+\sum_{i=1}^{n}\alpha_{n+i}\langle\mathcal{L}_{2}(1,0,0),\psi_{i}\rangle,\sum_{i=1}^{n}\alpha_{2n+i}\langle\mathcal{L}_{3}(1,0,0),\psi_{i}\rangle,\\\
\displaystyle\frac{\partial\mathfrak{L}}{\partial
b_{2}}&=0\rightarrow\displaystyle\sum_{i=1}^{n}\alpha_{i}\langle\mathcal{L}_{1}(0,1,0),\psi_{i}\rangle+\sum_{i=1}^{n}\alpha_{n+i}\langle\mathcal{L}_{2}(0,1,0),\psi_{i}\rangle,\sum_{i=1}^{n}\alpha_{2n+i}\langle\mathcal{L}_{3}(0,1,0),\psi_{i}\rangle,\\\
\displaystyle\frac{\partial\mathfrak{L}}{\partial
b_{3}}&=0\rightarrow\displaystyle\sum_{i=1}^{n}\alpha_{i}\langle\mathcal{L}_{1}(0,0,1),\psi_{i}\rangle+\sum_{i=1}^{n}\alpha_{n+i}\langle\mathcal{L}_{2}(0,0,1),\psi_{i}\rangle,\sum_{i=1}^{n}\alpha_{2n+i}\langle\mathcal{L}_{3}(0,0,1),\psi_{i}\rangle,\\\
\end{cases}$ $\displaystyle\frac{\partial\mathfrak{L}}{\partial\alpha_{k}}$
$\displaystyle=0\rightarrow\begin{cases}\displaystyle\sum_{j=1}^{d}w_{j}\langle\mathcal{L}_{1}(\varphi_{j},0,0)-f_{1},\psi_{k}\rangle+\sum_{j=1}^{d}w_{d+j}\langle\mathcal{L}_{1}(0,\varphi_{j},0)-f_{1},\psi_{k}\rangle+\\\
\displaystyle\sum_{j=1}^{d}w_{2d+j}\langle\mathcal{L}_{1}(0,0,\varphi_{j})-f_{1},\psi_{k}\rangle=e_{k}&k=1,2,\ldots,n.\\\
\displaystyle\sum_{j=1}^{d}w_{j}\langle\mathcal{L}_{2}(\varphi_{j},0,0)-f_{2},\psi_{k}\rangle+\sum_{j=1}^{d}w_{d+j}\langle\mathcal{L}_{2}(0,\varphi_{j},0)-f_{2},\psi_{k}\rangle+\\\
\displaystyle\sum_{j=1}^{d}w_{2d+j}\langle\mathcal{L}_{2}(0,0,\varphi_{j})-f_{2},\psi_{k}\rangle=e_{k}&k=n+1,n+2,\ldots,2n.\\\
\displaystyle\sum_{j=1}^{d}w_{j}\langle\mathcal{L}_{3}(\varphi_{j},0,0)-f_{3},\psi_{k}\rangle+\sum_{j=1}^{d}w_{d+j}\langle\mathcal{L}_{3}(0,\varphi_{j},0)-f_{3},\psi_{k}\rangle+\\\
\displaystyle\sum_{j=1}^{d}w_{2d+j}\langle\mathcal{L}_{3}(0,0,\varphi_{j})-f_{3},\psi_{k}\rangle=e_{k}&k=2n+1,2n+2,\ldots,3n.\\\
\end{cases}$

By defining

$\displaystyle A^{(1)}_{i,j,k}$
$\displaystyle=\langle\mathcal{L}_{1}(\varphi_{i},0,0),\psi_{j}\rangle,A^{(2)}_{i,j,k}=\langle\mathcal{L}_{2}(\varphi_{i},0,0),\psi_{j}\rangle,A^{(3)}_{i,j,k}=\langle\mathcal{L}_{3}(\varphi_{i},0,0),\psi_{j}\rangle,$
$\displaystyle A^{(4)}_{i,j,k}$
$\displaystyle=\langle\mathcal{L}_{1}(0,\varphi_{i},0),\psi_{j}\rangle,A^{(5)}_{i,j,k}=\langle\mathcal{L}_{2}(0,\varphi_{i},0),\psi_{j}\rangle,A^{(6)}_{i,j,k}=\langle\mathcal{L}_{3}(0,\varphi_{i},0),\psi_{j}\rangle,$
$\displaystyle A^{(7)}_{i,j,k}$
$\displaystyle=\langle\mathcal{L}_{1}(0,0,\varphi_{i}),\psi_{j}\rangle,A^{(8)}_{i,j,k}=\langle\mathcal{L}_{2}(0,0,\varphi_{i}),\psi_{j}\rangle,A^{(9)}_{i,j,k}=\langle\mathcal{L}_{3}(0,0,\varphi_{i}),\psi_{j}\rangle,$
$\displaystyle B^{(1)}_{i,j,k}$
$\displaystyle=\langle\mathcal{L}_{1}(1,0,0),\psi_{j}\rangle,B^{(2)}_{i,j,k}=\langle\mathcal{L}_{2}(1,0,0),\psi_{j}\rangle,B^{(3)}_{i,j,k}=\langle\mathcal{L}_{3}(1,0,0),\psi_{j}\rangle,$
$\displaystyle B^{(4)}_{i,j,k}$
$\displaystyle=\langle\mathcal{L}_{1}(0,1,0),\psi_{j}\rangle,B^{(5)}_{i,j,k}=\langle\mathcal{L}_{2}(0,1,0),\psi_{j}\rangle,B^{(6)}_{i,j,k}=\langle\mathcal{L}_{3}(0,1,0),\psi_{j}\rangle,$
$\displaystyle B^{(7)}_{i,j,k}$
$\displaystyle=\langle\mathcal{L}_{1}(0,0,1),\psi_{j}\rangle,B^{(8)}_{i,j,k}=\langle\mathcal{L}_{2}(0,0,1),\psi_{j}\rangle,B^{(9)}_{i,j,k}=\langle\mathcal{L}_{3}(0,0,1),\psi_{j}\rangle,$

and

$\displaystyle Z$ $\displaystyle=\begin{bmatrix}A^{(1)}&A^{(2)}&A^{(3)}\\\
A^{(4)}&A^{(5)}&A^{(6)}\\\ A^{(7)}&A^{(8)}&A^{(9)}\end{bmatrix},$
$\displaystyle V$ $\displaystyle=\begin{bmatrix}B^{(1)}&B^{(2)}&B^{(3)}\\\
B^{(4)}&B^{(5)}&B^{(6)}\\\ B^{(7)}&B^{(8)}&B^{(9)}\end{bmatrix},$

the relation (25) can be reformulated as:

$\displaystyle\begin{cases}Z\boldsymbol{\alpha}=\mathbf{w}\\\
\mathbf{e}=-\boldsymbol{\alpha}/\gamma\\\ b=V^{T}\boldsymbol{\alpha}\\\
Z^{T}\mathbf{w}-\mathbf{e}=y.\end{cases}$

Eliminating $\mathbf{w}$ and $\mathbf{e}$ yields

$\displaystyle\left[\begin{array}[]{c|c}\mathbf{0}&V^{T}\\\ \hline\cr
V&\Omega+I/\gamma\end{array}\right]\left[\begin{array}[]{c}b\\\
\hline\cr\boldsymbol{\alpha}\end{array}\right]=\left[\begin{array}[]{c}0\\\
\hline\cr y\end{array}\right],$ (26)

where

$\displaystyle\boldsymbol{\alpha}$
$\displaystyle=\left[\alpha_{1},\ldots,\alpha_{3n}\right]^{T},$ $\displaystyle
y$ $\displaystyle=\left[\langle f_{1},\psi_{1}\rangle,\langle
f_{1},\psi_{2}\rangle,\ldots,\langle f_{1},\psi_{n}\rangle,\langle
f_{2},\psi_{1}\rangle,\langle f_{2},\psi_{2}\rangle,\ldots,\langle
f_{2},\psi_{n}\rangle,\langle f_{3},\psi_{1}\rangle,\langle
f_{3},\psi_{2}\rangle,\ldots,\langle f_{3},\psi_{n}\rangle\right]^{T},$

and

$\Omega=Z^{T}Z=\begin{bmatrix}A^{T}&C^{T}\\\
B^{T}&D^{T}\end{bmatrix}\begin{bmatrix}A&B\\\ C&D\end{bmatrix}.$ (27)

The kernel trick also appears at each block of matrix $\Omega$. As a result,
the approximated solution in the dual form can be computed using:

$\displaystyle\tilde{u}_{1}(x)$
$\displaystyle=\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{1}(x,x_{i})+\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{2}(x,x_{i})+\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{3}(x,x_{i})+b_{1},$
$\displaystyle\tilde{u}_{2}(x)$
$\displaystyle=\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{4}(x,x_{i})+\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{5}(x,x_{i})\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{6}(x,x_{i})+b_{2},$
$\displaystyle\tilde{u}_{3}(x)$
$\displaystyle=\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{7}(x,x_{i})+\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{8}(x,x_{i})\sum_{i=1}^{n}\alpha_{i}\widetilde{K}_{9}(x,x_{i})+b_{3},$

where

$\displaystyle\widetilde{K}_{1}(x,x_{i})=\langle\langle\mathcal{L}_{1}(\varphi,0,0),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{1}(K(x,t),0,0),\psi_{i}\rangle,$
$\displaystyle\widetilde{K}_{2}(x,x_{i})=\langle\langle\mathcal{L}_{2}(\varphi,0,0),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{2}(K(x,t),0,0),\psi_{i}\rangle,$
$\displaystyle\widetilde{K}_{3}(x,x_{i})=\langle\langle\mathcal{L}_{3}(\varphi,0,0),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{3}(K(x,t),0,0),\psi_{i}\rangle,$
$\displaystyle\widetilde{K}_{4}(x,x_{i})=\langle\langle\mathcal{L}_{1}(0,\varphi,0),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{1}(0,K(x,t),0),\psi_{i}\rangle,$
$\displaystyle\widetilde{K}_{5}(x,x_{i})=\langle\langle\mathcal{L}_{2}(0,\varphi,0),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{2}(0,K(x,t),0),\psi_{i}\rangle,$
$\displaystyle\widetilde{K}_{6}(x,x_{i})=\langle\langle\mathcal{L}_{3}(0,\varphi,0),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{3}(0,K(x,t),0),\psi_{i}\rangle,$
$\displaystyle\widetilde{K}_{7}(x,x_{i})=\langle\langle\mathcal{L}_{1}(0,0,\varphi),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{1}(0,0,K(x,t)),\psi_{i}\rangle,$
$\displaystyle\widetilde{K}_{8}(x,x_{i})=\langle\langle\mathcal{L}_{2}(0,0,\varphi),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{2}(0,0,K(x,t)),\psi_{i}\rangle,$
$\displaystyle\widetilde{K}_{9}(x,x_{i})=\langle\langle\mathcal{L}_{3}(0,0,\varphi),\psi_{i}\rangle,\varphi\rangle=\langle\mathcal{L}_{3}(0,0,K(x,t)),\psi_{i}\rangle.$

Hence, the proposed approach formulates the solutions (22) using kernel
functions $(K(x,t))$ which are replaced by Legendre polynomials [16]. Also, it
would be possible to employ $\boldsymbol{\alpha}$ to calculate $\mathbf{w}$,
and by utilizing Formula (22), one could approximate the solution.

######  Remark 1.

Although the choice of basis functions can differ for each approximated
function, for simplicity, uniform functions are employed. In this proposed
methodology, Legendre polynomials (15) of degree one, denoted as $P_{1}(x)=x$,
are utilized. Additionally, considering that $P_{0}(x)=1$, the variable $b$
(bias term) is omitted from this formulation.

######  Remark 2.

Since the Legendre function is defined in the interval $(-1,1)$ in the formula
(14), we built the kernel of the network by the shifted Legendre polynomials
using (16).

######  Remark 3.

The collocation method considers the test function as
$\psi_{j}(x)=\delta(x-x_{j})$, where $\delta$ is the Dirac delta function.
According to the characteristic of Dirac’s Delta function, in the formula (20)
we have
$\langle\mathcal{N}_{i}(\tilde{u}_{1},\tilde{u}_{2},\ldots,\tilde{u}_{k})-f_{i},\delta(x-x_{j})\rangle=\mathcal{N}_{i}(\tilde{u}_{1}(x_{j}),\tilde{u}_{2}(x_{j}),\ldots,\tilde{u}_{k}(x_{j}))-f_{i}(x_{j})$.

This methodology is termed Collocation Least-Squares Support Vector Regression
(CLS-SVR). In this algorithm, we employ the roots of Legendre polynomials of
degree $m$ as the training data for the CLS-SVR learning process.

##  4 Numerical results

To show the efficiency of the proposed method, this section presents a
numerical solution to various well-known differential-algebraic equations.
This example includes DAEs, VI-DAE of fractional order, linear and nonlinear
FDAEs, and a system of partial differential-algebraic equations. All the
results have been obtained by running the algorithm using Maple 20. The error
criterion of the presented tables is the relative error defined by:

$\mathbb{E}_{u}=|\frac{u(t)-\tilde{u}(t)}{u(t)}|,$ (28)

where $u(t)$ is the exact solution and $\tilde{u}(t)$ is the predicted answer
of the algorithm. In all these examples, the value of the regularization
coefficient $\gamma$ is chosen between $10$ and $10^{3}$. These differences
are caused by the structure of the problem and the limitations of Maple. In
addition to the numerical tables and relative error for each example, the
graph of the absolute error is given.

######  Example 1.

As the first example, consider the following differential-algebraic equation

(a)

(b)

(c)

Figure 1: (a) $|u_{1}(t)-\tilde{u}_{1}(t)|$, (b)
$|u_{2}(t)-\tilde{u}_{2}(t)|$, and (c) $|u_{3}(t)-\tilde{u}_{3}(t)|$ for
example 1 with $m=4$. Table 2: Exact, approximate solution, and relative error
for $m=10$ on $\left[0,1\right]$ for example 1

.  $t$ $u_{1}(t)$ $\tilde{u}_{1}(t)$ $\mathbb{E}_{u_{1}(t)}$ 0.2 0.0397338
0.0397337 $1.8\times 10^{-6}$ 0.4 0.15576733 0.15576736 $1.6\times 10^{-7}$
0.6 0.3387854 0.3387855 $1.2\times 10^{-7}$ 0.8 0.5738848 0.5738847 $1.7\times
10^{-7}$ 1 0.84147098 0.84147097 $1.6\times 10^{-8}$ $u_{2}(t)$
$\tilde{u}_{2}(t)$ $\mathbb{E}_{u_{2}(t)}$ 0.20271 0.20270 $3.5\times 10^{-7}$
0.422793 0.422792 $1.6\times 10^{-6}$ 0.684136 0.684137 $1.5\times 10^{-6}$
1.02963855 1.02963859 $3.6\times 10^{-8}$ 1.5574077 1.5574078 $8.9\times
10^{-8}$ $u_{3}(t)$ $\tilde{u}_{3}(t)$ $\mathbb{E}_{u_{3}(t)}$ 0.196013
0.196017 $2.3\times 10^{-5}$ 0.3684243 0.3684244 $1.2\times 10^{-7}$ 0.495201
0.495196 $9.8\times 10^{-6}$ 0.55736 0.55737 $1.2\times 10^{-5}$ 0.540302
0.540293 $1.7\times 10^{-5}$

$\begin{cases}u^{\prime}_{1}(t)=u_{1}(t)-u_{3}(t)u_{2}(t)+sin(t)+tcos(t),\\\
u^{\prime}_{2}(t)=tu_{3}(t)+u^{2}_{1}(t)+sec^{2}(t)-t^{2}(cos(t)+sin^{2}(t)),\\\
0=u_{1}(t)-u_{3}(t)+t(cos(t)-sin(t)),\quad t\in[0,1]\\\ \end{cases}$ (29)

with initial conditions $u_{1}(0)=u_{2}(0)=u_{3}(0)=0$ and the exact solutions
as $u_{1}(t)=tsin(t),u_{2}(t)=tan(t)$ and $u_{3}(t)=tcos(t)$. For this case,
we set $m=10$ and $\gamma=10^{2}$ over the interval $\left[0,1\right]$. The
relative error value obtained via the CLS-SVR method for each $u_{i}$ is
denoted as $\mathbb{E}_{u_{i}}$. In Table 2, the error of this case is
presented. In Figure 1, the absolute error of the $u_{1}(t)$, $u_{2}(t)$, and
$u_{3}(t)$ is plotted.

######  Example 2.

In this example, we consider a fractional order integro-differential algebraic
equation :

$\begin{cases}D^{\frac{1}{2}}u_{1}(t)=+t\int_{0}^{t}u_{1}(s)ds+(1+t)\int_{0}^{t}u_{2}(s)ds+(\frac{3t\sqrt{\pi}}{4}-\frac{2t^{\frac{7}{2}}}{5}-\frac{2(1+t)t^{\frac{5}{2}}}{5}),\\\
0=\int_{0}^{t}(1+s)u_{1}(s)ds+\int_{0}^{t}u_{2}(s)ds+(\frac{2t^{\frac{5}{2}}(5t+7)}{35}+\frac{2t^{\frac{5}{2}}}{5}),\quad
t\in[0,1],\end{cases}$ (30)

with zero initial conditions and the exact solutions
$u_{1}(t)=u_{2}(t)=t\sqrt{t}$.

Table 3: Exact, approximate solution, and relative error for $m=14$ on
$\left[0,1\right]$ for example 2

.  $t$ $u_{1}(t)$ $\tilde{u}_{1}(t)$ $\mathbb{E}_{u_{1}(t)}$ 0.2 0.089442
0.089452 $1.1\times 10^{-4}$ 0.4 0.252982 0.252989 $2.8\times 10^{-5}$ 0.6
0.464758 0.464760 $5.5\times 10^{-6}$ 0.8 0.715541 0.715542 $1.5\times
10^{-6}$ 1 1 1.0000007 $7.3\times 10^{-7}$ $u_{2}(t)$ $\tilde{u}_{2}(t)$
$\mathbb{E}_{u_{2}(t)}$ 0.089442 0.089690 $2.7\times 10^{-3}$ 0.252982
0.253186 $8.0\times 10^{-4}$ 0.464758 0.465094 $7.3\times 10^{-4}$ 0.715541
0.715935 $5.5\times 10^{-4}$ 1 1.01319 $1.3\times 10^{-2}$

(a)

(b)

Figure 2: (a) $|u_{1}(t)-\tilde{u}_{1}(t)|$ and (b)
$|u_{2}(t)-\tilde{u}_{2}(t)|$ for example 2 with $m=14$.

The provided equations belong to the category of Volterra integro-differential
algebraic equations of fractional order. This study employs the Caputo
derivative (10) to account for fractional derivatives in the differential
equations. To address this equation, the proposed model is trained using
parameters $m=14$ and $\gamma=10^{3}$ across the interval $[0,1]$. Details of
the error analysis resulting from the program’s execution are provided in
Table 3 and Figure 2.

######  Example 3.

In the third problem, we calculate the solution of a linear fractional
differential-algebraic equation

$\begin{cases}D^{\frac{1}{2}}u_{1}(t)+2u_{1}(t)-\frac{\Gamma(\frac{7}{2})}{\Gamma(3)}u_{2}(t)+u_{3}(t)=2t^{\frac{5}{2}}+sin(t),\\\
D^{\frac{1}{2}}u_{2}(t)+u_{2}(t)+u_{3}(t)=\frac{\Gamma(3)}{\Gamma(\frac{5}{2})}t^{\frac{3}{2}}+t^{2}+sin(t),\\\
2u_{1}(t)+u_{2}(t)-u_{3}(t)=2t^{\frac{5}{2}}+t^{2}-sin(t),\quad t\in[0,1]\\\
\end{cases}$ (31)

in the context of zero initial conditions. The corresponding exact solutions
are given by $u_{1}(t)=t^{\frac{5}{2}}$, $u_{2}(t)=t^{2}$, and
$u_{3}(t)=\sin(t)$.

(a)

(b)

(c)

Figure 3: (a) $|u_{1}(t)-\tilde{u}_{1}(t)|$, (b)
$|u_{2}(t)-\tilde{u}_{2}(t)|$, and (c) $|u_{3}(t)-\tilde{u}_{3}(t)|$ for
example 3 with $m=10$. Table 4: Exact, approximate solution, and relative
error for $m=10$ on $\left[0,1\right]$ for example 3

.  $t$ $u_{1}(t)$ $\tilde{u}_{1}(t)$ $\mathbb{E}_{u_{1}(t)}$ 0.2 0.0178885
0.0178883 $1.3\times 10^{-5}$ 0.4 0.1011928 0.10119127 $1.5\times 10^{-5}$ 0.6
0.2788548 0.27885540 $2.1\times 10^{-6}$ 0.8 0.57243340 0.5724344 $1.7\times
10^{-6}$ 1 1 0.9999 $1.0\times 10^{-6}$ $u_{2}(t)$ $\tilde{u}_{2}(t)$
$\mathbb{E}_{u_{2}(t)}$ 0.04 0.0400003 $8.2\times 10^{-6}$ 0.16 0.159999
$7.4\times 10^{-7}$ 0.36 0.359999 $9.9\times 10^{-8}$ 0.64 0.6400002
$3.8\times 10^{-7}$ 1 0.99999968 $3.1\times 10^{-7}$ $u_{3}(t)$
$\tilde{u}_{3}(t)$ $\mathbb{E}_{u_{3}(t)}$ 0.198669 0.198667 $8.8\times
10^{-6}$ 0.389418 0.389419 $2.1\times 10^{-6}$ 0.564642 0.564643 $1.4\times
10^{-6}$ 0.717356 0.717354 $2.1\times 10^{-6}$ 0.841470 0.841477 $8.4\times
10^{-6}$

Figure 3 illustrates the difference between the estimated solution and the
exact solution for each $u_{i}(t)$. The relative error is computed across
various intervals within the range of $[0,1]$, as indicated in Table 4. We
compare our method with the waveform relaxation scheme proposed in [8]. In
Table 5, we present a comparison based on the $\ell_{2}$ norm error and the
highest iteration of $N=20$, using the method described in [8].

Table 5: Comparison of the obtained errors between our method and the method [8] for Example 3 | LS-SVR method |  method [8]
---|---|---
$\|u_{1}(t)-\tilde{u}_{1}(t)\|_{2}$ | $2.2\times 10^{-6}$ | $\approx 10^{-3}$
$\|u_{2}(t)-\tilde{u}_{2}(t)\|_{2}$ | $5.3\times 10^{-7}$ | $\approx 10^{-3}$
$\|u_{3}(t)-\tilde{u}_{3}(t)\|_{2}$ | $7.4\times 10^{-6}$ | $\approx 10^{-3}$

######  Example 4.

Consider the following nonlinear fractional differential algebraic equations

$\begin{cases}D^{\frac{1}{2}}u_{1}(t)+u_{1}(t)u_{2}(t)-u_{3}(t)=\frac{\Gamma(4)}{\Gamma(\frac{7}{2})}t^{\frac{5}{2}}+2t^{4}+t^{7}-e^{t}-tsin(t),\\\
D^{\frac{1}{2}}u_{2}(t)-\frac{\Gamma(5)}{\Gamma(\frac{9}{2})}t^{\frac{1}{2}}u_{1}(t)+2u_{2}(t)+u_{1}(t)u_{3}(t)={\frac{2}{\Gamma(\frac{3}{2})}}t^{\frac{1}{2}}+4t+2t^{4}+t^{3}e^{t}+t^{4}sin(t),\\\
u^{2}_{1}(t)-t^{2}u_{2}(t)+u_{3}(t)=e^{t}+2t^{3}+tsin(t),\quad t\in[0,1],\\\
\end{cases}$ (32)

with the consistent initial conditions $u_{1}(0)=0,u_{2}(0)=0,$ and
$u_{3}(0)=1$. The equation has analytical solutions
$u_{1}(t)=t^{3},u_{2}(t)=2t+t^{4},u_{3}(t)=e^{t}+tsin(t)$. Table 7 calculates
the numerical results and the equation’s absolute error value. Comparison of
our reported results with those obtained in Table 6 confirms the superiority
and reliability of our proposed scheme over the method presented in [8, 10]
with the highest iteration of $N=16$.

Table 6: Comparison of the obtained errors between our method and methods in [8, 10] for Example 4 | LS-SVR method |  method [8] |  method [10]
---|---|---|---
$\|u_{1}(t)-\tilde{u}_{1}(t)\|_{2}$ | $2.4\times 10^{-14}$ | $\approx 10^{-3}$ | $2.3\times 10^{-7}$
$\|u_{2}(t)-\tilde{u}_{2}(t)\|_{2}$ | $6.3\times 10^{-15}$ | $\approx 10^{-3}$ | $9.0\times 10^{-5}$
$\|u_{3}(t)-\tilde{u}_{3}(t)\|_{2}$ | $9.8\times 10^{-14}$ | $\approx 10^{-4}$ | $4.6\times 10^{-5}$

(a)

(b)

(c)

Figure 4: (a) $|u_{1}(t)-\tilde{u}_{1}(t)|$, (b)
$|u_{2}(t)-\tilde{u}_{2}(t)|$, and (c) $|u_{3}(t)-\tilde{u}_{3}(t)|$ for
example 4 with $m=10$. Table 7: Exact, approximate solution, and relative
error for $m=10$ on $\left[0,1\right]$ for example 4

.  $t$ $u_{1}(t)$ $\tilde{u}_{1}(t)$ $\mathbb{E}_{u_{1}(t)}$ 0.2 0.008
0.007999 $2.8\times 10^{-13}$ 0.4 0.064 0.063999 $2.7\times 10^{-13}$ 0.6
0.216 0.216000 $3.3\times 10^{-14}$ 0.8 0.512 0.512000 $2.8\times 10^{-14}$ 1
1 0.999999 $1.2\times 10^{-15}$ $u_{2}(t)$ $\tilde{u}_{2}(t)$
$\mathbb{E}_{u_{2}(t)}$ 0.4016 0.4016000 $1.9\times 10^{-15}$ 0.8256 0.8256000
$7.6\times 10^{-16}$ 1.3296 1.3295999 $1.1\times 10^{-15}$ 2.0096 2.0095999
$3.0\times 10^{-15}$ 3 2.9999999 $6.4\times 10^{-17}$ $u_{3}(t)$
$\tilde{u}_{3}(t)$ $\mathbb{E}_{u_{3}(t)}$ 1.2611366243191 1.2611366243192
$1.5\times 10^{-14}$ 1.64759203456473 1.64759203456476 $2.0\times 10^{-14}$
2.16090428442753 2.16090428442750 $1.0\times 10^{-14}$ 2.79942580121208
2.79942580121209 $2.9\times 10^{-15}$ 3.559752813266 3.559752813267 $2.4\times
10^{-14}$

######  Example 5.

Consider the partial differential-algebraic equation

$\begin{cases}D_{t}u_{2}(x,t)+D_{t}u_{3}(x,t)=-\frac{1}{2}x^{2}e^{-\frac{1}{2}t}+x^{2}cos(t),\\\
2D_{t}u_{1}(x,t)-D_{t}u_{2}(x,t)-D_{t}u_{3}(x,t)-u_{2}(x,t)=-2x^{2}e^{-t}-\frac{x^{2}e^{\frac{-t}{2}}}{2}-x^{2}cos(t),\\\
-D_{xx}u_{3}(x,t)+u_{3}(x,t)=-2sin(t)+x^{2}sin(t),\quad x\in[-0.5,0.5],\quad
t\in[0,1],\\\ \end{cases}$ (33)

with the initial conditions:

$u_{1}(x,0)=x^{2},D_{t}u_{1}(x,0)=-x^{2},u_{2}(x,0)=x^{2},D_{t}u_{2}(x,0)=\frac{-x^{2}}{2},u_{3}(x,0)=0,D_{t}u_{3}(x,0)=x^{2},$

and whose analytical solution is $u(x,t)=\begin{bmatrix}x^{2}e^{-t}\\\
x^{2}e^{\frac{-t}{2}}\\\ x^{2}sin(t)\end{bmatrix}$. In Table 8 the relative
error is calculated for a number of points in the interval $[0,1]$. the
absolute error with $m=6$ and $\gamma=100$ is in Figure 5. Also, the relative
error for various $m$ has been calculated in Table 9.

(a)

(b)

(c)

Figure 5: (a) $|u_{1}(x,t)-\tilde{u}_{1}(x,t)|$, (b)
$|u_{2}(x,t)-\tilde{u}_{2}(x,t)|$, and (c) $|u_{3}(x,t)-\tilde{u}_{3}(x,t)|$
for example 5 with m=6. Table 8: Exact, approximate solution, and relative
error of example 5 with $m=6$.

$t$ $x$ $u_{1}(x,t)$ $\tilde{u}_{1}(x,t)$ $\mathbb{E}_{u_{1}(x,t)}$ 0.02 0.02 0.0003920794 0.0003920795 $3.2\times 10^{-7}$ 0.04 0.04 0.0015372631 0.0015372635 $2.9\times 10^{-7}$ 0.06 0.06 0.00339035232 0.0033903531 $2.5\times 10^{-7}$ 0.08 0.08 0.005907944 0.005907945 $2.1\times 10^{-7}$ 0.1 0.1 0.009048374 0.009048375 $1.8\times 10^{-7}$ |  $u_{2}(x,t)$ $\tilde{u}_{2}(x,t)$ $\mathbb{E}_{u_{2}(x,t)}$ 0.00039601 0.00039602 $7.8\times 10^{-6}$ 0.00156831 0.00156832 $6.8\times 10^{-6}$ 0.00349360 0.00349362 $5.8\times 10^{-6}$ 0.00614905 0.00614908 $4.9\times 10^{-6}$ 0.00951229 0.0095123 $4.2\times 10^{-6}$ |  $u_{3}(x,t)$ $\tilde{u}_{3}(x,t)$ $\mathbb{E}_{u_{3}(x,t)}$ 0.0000079994 0.000007996 $3.9\times 10^{-4}$ 0.00006398 0.00006397 $1.6\times 10^{-4}$ 0.00021587 0.00021584 $9.5\times 10^{-5}$ 0.00051145 0.00051142 $6.0\times 10^{-5}$ 0.0009983 0.0009982 $4.0\times 10^{-5}$
---|---|---

Table 9: Relative error of example 5 with different $m$ for $(x,t)=(.02,0.02),(.04,.04),(.06,.06),(.08,.08),and(.1,.1)$ $m$ | $\mathbb{E}_{u_{1}(x,t)}$ | $\mathbb{E}_{u_{2}(x,t)}$ | $\mathbb{E}_{u_{3}(x,t)}$
---|---|---|---
$m=8$ |  $2.9\times 10^{-9}$ $2.1\times 10^{-9}$ $1.4\times 10^{-9}$ $8.7\times 10^{-10}$ $4.0\times 10^{-10}$ |  $3.9\times 10^{-8}$ $1.9\times 10^{-8}$ $5.1\times 10^{-9}$ $5.2\times 10^{-9}$ $1.2\times 10^{-8}$ |  $1.9\times 10^{-6}$ $4.8\times 10^{-7}$ $8.3\times 10^{-8}$ $6.3\times 10^{-8}$ $1.1\times 10^{-7}$
$m=10$ |  $5.2\times 10^{-13}$ $3.6\times 10^{-13}$ $2.5\times 10^{-13}$ $1.7\times 10^{-13}$ $1.1\times 10^{-13}$ |  $2.8\times 10^{-11}$ $1.7\times 10^{-11}$ $9.5\times 10^{-12}$ $4.6\times 10^{-12}$ $1.5\times 10^{-12}$ |  $1.3\times 10^{-9}$ $4.2\times 10^{-10}$ $1.5\times 10^{-10}$ $5.5\times 10^{-11}$ $1.5\times 10^{-11}$
$m=12$ |  $1.7\times 10^{-15}$ $1.0\times 10^{-15}$ $5.7\times 10^{-16}$ $3.0\times 10^{-16}$ $1.5\times 10^{-16}$ |  $1.2\times 10^{-13}$ $5.5\times 10^{-14}$ $1.8\times 10^{-14}$ $1.4\times 10^{-15}$ $5.0\times 10^{-15}$ |  $6.1\times 10^{-12}$ $1.3\times 10^{-12}$ $3.0\times 10^{-13}$ $1.6\times 10^{-14}$ $4.8\times 10^{-14}$

##  5 Conclusion

The primary aim of this research is to establish a cohesive relationship among
a machine learning model, weighted residual methods, and Legendre’s orthogonal
polynomials to effectively address DAEs. We adopted a system of differential-
algebraic equations and employed LS-SVR to approximate their solutions. We
represented the solution as a linear combination of orthogonal Legendre
polynomials with unknown coefficients and subsequently formulated the residual
function. This residual function was aligned with the principles of weighted
residual methods by establishing its inner product with the Dirac delta
function. The resulting system of equations was then adopted as the
constraints for the quadratic programming problem. This approach presents
notable advantages in contrast to conventional methods, including the
utilization of regularization techniques for ill-conditioned problems,
attainment of a singular solution with exponential convergence rate, and
formulation of a quadratic programming problem leading to a symmetric positive
definite matrix within a dual space. These attributes collectively contribute
to the inherent simplicity of this technique, rendering it suitable for
solving diverse categories of DAEs. Through the examples presented, we have
applied LS-SVR for solving the Volterra integral DAE of fractional order,
partial DAE, and fractional DAEs with high accuracy for solving equations.
Nonetheless, this approach encounters challenges in solving DAEs with non-
periodic boundary conditions. Addressing the issue within the current
framework is complex and requires additional research to improve and broaden
the method. Within this model, we achieved precise outcomes through the
application of widely recognized Legendre orthogonal functions as the model’s
kernel. For future investigations, alternative kernel functions like
Chebyshev, Hermite, and Gegenbauer polynomials hold potential for exploration.
Furthermore, the scope of this method can be expanded to encompass the
resolution of other differential equation categories.

##  6 Declarations

Declarations of interest: none

##  7 Data availability statement

not applicable.

## References

  * Aghaei and Parand [2023] A. A. Aghaei and K. Parand.  Hyperparameter optimization of orthogonal functions in the numerical solution of differential equations.  _arXiv preprint arXiv:2304.14088_ , 2023.
  * Ahmed [2020] H. F. Ahmed.  Numerical study on factional differential-algebraic systems by means of Chebyshev Pseudo spectral method.  _Journal of Taibah University for Science_ , 14(1):1023–1032, 2020.
  * Benhammouda and Vazquez-Leal [2014] B. Benhammouda and H. Vazquez-Leal.  Analytical solutions for systems of partial differential–algebraic equations.  _SpringerPlus_ , 3(1):1–9, 2014.
  * Butcher and Chan [1998] J. Butcher and R. Chan.  Efficient Runge-Kutta integrators for index-2 differential algebraic equations.  _Mathematics of Computation_ , 67(223):1001–1021, 1998.
  * Caputo [1967] M. Caputo.  Linear models of dissipation whose Q is almost frequency independent—II.  _Geophysical Journal International_ , 13(5):529–539, 1967.
  * Çelik and Bayram [2003] E. Çelik and M. Bayram.  On the numerical solution of differential-algebraic equations by Padé series.  _Applied mathematics and computation_ , 137(1):151–160, 2003.
  * Celik et al. [2006] E. Celik, M. Bayram, and T. Yeloglu.  Solution of differential-algebraic equations (DAEs) by adomian decomposition method.  _International Journal Pure & Applied Mathematical Sciences_, 3(1):93–100, 2006.
  * Ding and Jiang [2014] X.-L. Ding and Y.-L. Jiang.  Waveform relaxation method for fractional differential-algebraic equations.  _Fractional Calculus and Applied Analysis_ , 17:585–604, 2014.
  * Foupouagnigni and Koepf [2020] M. Foupouagnigni and W. Koepf.  _Orthogonal Polynomials: 2nd AIMS-Volkswagen Stiftung Workshop, Douala, Cameroon, 5-12 October, 2018_.  Springer Nature, 2020.
  * Ghanbari et al. [2018] F. Ghanbari, K. Ghanbari, and P. Mokhtary.  Generalized Jacobi–Galerkin method for nonlinear fractional differential algebraic equations.  _Computational and Applied Mathematics_ , 37:5456–5475, 2018.
  * Günerhan and Çelik [2020] H. Günerhan and E. Çelik.  Analytical and approximate solutions of fractional partial differential-algebraic equations.  _Applied Mathematics and Nonlinear Sciences_ , 5(1):109–120, 2020.
  * Hajimohammadi et al. [2020] Z. Hajimohammadi, F. Baharifard, and K. Parand.  A new numerical learning approach to solve general Falkner–Skan model.  _Engineering with Computers_ , pages 1–17, 2020.
  * Hosseini [2006] M. M. Hosseini.  Adomian decomposition method for solution of nonlinear differential algebraic equations.  _Applied mathematics and computation_ , 181(2):1737–1744, 2006.
  * Lu et al. [2020] Y. Lu, Q. Yin, H. Li, H. Sun, Y. Yang, and M. Hou.  Solving higher order nonlinear ordinary differential equations with least squares support vector machines.  _Journal of Industrial and Management Optimization_ , 16(3):1481–1502, 2020.
  * Miller and Ross [1993] K. S. Miller and B. Ross.  _An introduction to the fractional calculus and fractional differential equations_.  Wiley, 1993.
  * Parand et al. [2021] K. Parand, A. A. Aghaei, M. Jani, and A. Ghodsi.  A new approach to the numerical solution of Fredholm integral equations using least squares-support vector regression.  _Mathematics and Computers in Simulation_ , 180:114–128, 2021.
  * Shen et al. [2011] J. Shen, T. Tang, and L.-L. Wang.  _Spectral methods: algorithms, analysis and applications_ , volume 41.  Springer Science & Business Media, 2011.
  * Soltanian et al. [2009] F. Soltanian, S.-M. Karbassi, and M. Hosseini.  Application of He’s variational iteration method for solution of differential-algebraic equations.  _Chaos, Solitons & Fractals_, 41(1):436–445, 2009.
  * Soltanian et al. [2010] F. Soltanian, M. Dehghan, and S.-M. Karbassi.  Solution of the differential algebraic equations via homotopy perturbation method and their engineering applications.  _International Journal of Computer Mathematics_ , 87(9):1950–1974, 2010.
  * Suykens et al. [2002] J. A. Suykens, T. Van Gestel, J. De Brabanter, B. De Moor, and J. Vandewalle.  Least squares support vector machines, 2002.
  * Taheri et al. [2023] T. Taheri, A. A. Aghaei, and K. Parand.  Bridging machine learning and weighted residual methods for delay differential equations of fractional order.  _Applied Soft Computing_ , 149:110936, 2023.
  * Zolfaghari et al. [2021] R. Zolfaghari, J. Taylor, and R. J. Spiteri.  Structural analysis of integro-differential–algebraic equations.  _Journal of Computational and Applied Mathematics_ , 394:113568, 2021.
  * Zurigat et al. [2010] M. Zurigat, S. Momani, and A. Alawneh.  Analytical approximate solutions of systems of fractional algebraic–differential equations by homotopy analysis method.  _Computers & Mathematics with Applications_, 59(3):1227–1235, 2010.