# Physics-Informed Neural Networks for Discovering Localised Eigenstates in
Disordered Media

Liam Harcombe<EMAIL_ADDRESS>Quanling Deng
<EMAIL_ADDRESS>School of Computing, Australian National University,
Canberra, ACT 2601, Australia

###### Abstract

The Schrödinger equation with random potentials is a fundamental model for
understanding the behavior of particles in disordered systems. Disordered
media are characterised by complex potentials that lead to the localisation of
wavefunctions, also called Anderson localisation. These wavefunctions may have
similar scales of eigenenergies which poses difficulty in their discovery. It
has been a longstanding challenge due to the high computational cost and
complexity of solving the Schrödinger equation. Recently, machine-learning
tools have been adopted to tackle these challenges. In this paper, based upon
recent advances in machine learning, we present a novel approach for
discovering localised eigenstates in disordered media using physics-informed
neural networks (PINNs). We focus on the spectral approximation of
Hamiltonians in one dimension with potentials that are randomly generated
according to the Bernoulli, normal, and uniform distributions. We introduce a
novel feature to the loss function that exploits known physical phenomena
occurring in these regions to scan across the domain and successfully discover
these eigenstates, regardless of the similarity of their eigenenergies. We
present various examples to demonstrate the performance of the proposed
approach and compare it with isogeometric analysis.

###### keywords:

Schrödinger equation , Hamiltonian , Anderson localisation , eigenvalues and
eigenstates , neural networks , isogeometric analysis

††journal: arXiv

##  1 Introduction and Background

Partial differential equations (PDEs) are powerful tools for studying dynamic
systems and have applications across various fields. Despite their usefulness,
these equations are notoriously difficult to solve explicitly, leading to an
increased interest in numerical approximations. In addition to classical
numerical methods, recent research has shown the potential of using neural
networks to estimate solutions to differential equations [26, 27, 31]. There
are several advantages of using neural networks to solve PDEs over traditional
numerical methods.  Firstly, numerical errors are usually not accumulated
[33]. Secondly, they are more robust against the “curse of dimensionality,”
which refers to the difficulty of solving PDEs in high-dimensional spaces
using traditional numerical methods [19]. Furthermore, they can be trained to
solve a wide range of PDEs, including those with complex boundary conditions
or nonlinear operators [16].  Once trained, neural networks can quickly
produce solutions for new PDEs that are similar to those encountered during
training, potentially without requiring further adjustments or recalibration
of parameters [13]. Overall, neural networks can provide a powerful and
flexible tool for solving PDEs, especially in cases where traditional
numerical methods may be insufficient or computationally expensive.

Data-driven supervised networks [15] and data-free unsupervised networks [23]
have been shown to produce efficient approximations of differential eigenvalue
problems. Supervised networks aim to discover patterns in labelled datasets to
infer a general model that is able to predict examples outside of the given
dataset. The training process involves calculating approximations for given
inputs and computing how much these estimations deviate from the true outputs
using a loss function. Then, the networks work to minimise this loss function
by adjusting its estimation process. Artificial Neural Networks (ANNs) are
comprised of layers of interconnected nodes with weights and biases that
calculate the estimations. Convolutional Neural Networks (CNNs) are a newer
form of neural networks that are especially known for their ability to
identify patterns faster and more reliably than ANNs in solving eigenvalue
problems using labelled datasets [15]. They introduce a convolutional layer
built up of windows that slide around the input data, searching for patterns
in whole regions at a time. They are commonly used to pass information on to
ANNs and are optimised similarly. In situations where training data is
unavailable or difficult to obtain, unsupervised networks are employed as a
solution, operating without a labelled dataset. In these situations, more
attention is drawn to the design of the network’s loss function, as there is
no labelled data to compare predictions with [7].

In the case of finding the eigenstates of Hamiltonians, we leverage known
physical phenomena to design our loss function. This kind of model is an
example of a Physics Informed Neural Network (PINN, c.f., [35, 4]), one that
embeds the knowledge of physical laws governing the system into the learning
process to drive the network to an admissible approximation.

Differential eigenvalue problems occur in a wide range of problems in physics
and applied mathematics, such as quantum energy problems. In early works,
Lagaris et al. [27] proposed an ANN to discover the solutions to differential
eigenvalue problems, tested on various problems in quantum mechanics. More
recently, Finol et al. [15] presented a supervised ANN to solve eigenvalue
problems in mechanics, showing that CNNs were outperforming traditional ANNs
in contexts with labelled datasets. Sirignano and Spiliopoulos [36] developed
a deep learning network to accurately solve partial differential equations in
dimensions as high as 200. They also proposed a mesh-free algorithm, which was
desirable as meshes become infeasible in such high dimensions. Chen et al. [6]
adopted the emerging PINNs to solve representative inverse scattering problems
in photonic metamaterials and nano-optic technologies. Their method was also
mesh-free and was tested against numerical simulations based on the finite
element method. These ANN models are excellent at learning time series data.
The work [32] demonstrated its application in wave packet displacements in
localised and delocalised quantum states. Another work is [25], where they
adopted CNNs to solve a classification problem. The authors constructed a
supervised model learning from experimental data to distinguish between the
dynamics of an Anderson insulator and a many-body localised phase. Yang et al.
[38] proposed a combination of PINNs with Bayesian Neural Networks, which
solved both forward and inverse nonlinear problems, obtaining more accurate
predictions than PINNs in systems with large noises.

Jin et al. [23] have shown that PINNs can solve single and multiple square
well problems, calculating the eigenfunction and eigenvalue simultaneously.
Their network employs a scanning mechanism that pushes the network to search
for higher eigenvalues as the training process evolves, while the network
updates the eigenfunction prediction accordingly. They store found
eigenfunctions and exploit the orthogonality of wave functions to search for
higher eigenstates. Grubišić et al. [18] applied PINNs to identify localised
eigenstates of operators with random potential. They studied the effective
potential of the operator, whose local minima correspond to the different
eigenstates. They first built a PINN to solve these eigenvalue problems in one
dimension, then generalised to a deeper model that solves higher dimensional
problems. Effective potentials provide a neat way of locating the localised
eigenstates, however, are limited in the number of locations they can
identify, and have difficulty differentiating between eigenstates of identical
eigenvalue. These limitations occur in Bernoulli distributed potentials, which
we aim to solve with our model.

We build on the design of [23] to approximate these eigenstates for
Hamiltonians whose potential is distributed randomly, leading to localisation
of the solution to specific regions of the domain. In particular, the core of
our work is in finding these eigenstates where the eigenenergies are nearly
identical, a task that current models are unable to accomplish. This is most
prevalent with potentials distributed according to the Bernoulli distribution
but can occur in other distributions.

The rest of the paper is organised as follows. In Section 2, we state the
Schrödinger differential eigenvalue problem and present the isogeometric
analysis, followed by a discussion on the challenges of its spectral
approximation. In section 3, we propose our novel loss function design for the
PINN for solving the Schrödinger equation. Section 4 collects and discusses
various numerical tests to demonstrate the performance of the proposed method.
Concluding remarks are presented in section 5.

##  2 Problem Statement and Isogeometric Analysis

In this section, we first introduce the modeling problem, followed by a
presentation of a classic numerical method, namely isogeometric analysis, used
to solve the problem. We then provide a general discussion on the challenges
associated with the numerical computation of the model.

###  2.1 Problem Statement

We study the time-independent Schrödinger equation: Find the eigenpair $(E,u)$
such that

$\displaystyle\left[-\frac{\hbar^{2}}{2m}\Delta+V\right]u$
$\displaystyle=Eu\quad\text{in}\quad\Omega,$ (2.1) $\displaystyle u$
$\displaystyle=0\quad\text{on}\quad\partial\Omega,$

where $\Delta=\nabla^{2}$ is the Laplacian, $V=V(x)\in L^{2}(\Omega)$
specifies the potential and is a non-negative function, and
$\Omega\subset\mathbb{R}^{d},d=1,2,3,$ is a bounded open domain with Lipschitz
boundary $\partial\Omega$. Throughout the paper, we will focus on the case
with $d=1$, followed by an example of an extension to $d=2$ in Section 4.3.
The differential operator is referred to as the Hamiltonian, i.e.,
$\mathcal{H}=-\frac{\hbar^{2}}{2m}\Delta+V$. Here, $\hbar$ is the reduced
Planck constant and $m$ is the mass of the particle. $\frac{\hbar^{2}}{2m}$ is
a constant that we can divide throughout equation 2.1 and absorb into $V$ and
$E$. Thus, finding the eigenpairs in equation 2.1 is equivalent to finding
those for the following equation, up to a scalar:

$\displaystyle-\Delta u+Vu$ $\displaystyle=Eu\quad\text{in}\quad\Omega,$ (2.2)
$\displaystyle u$ $\displaystyle=0\quad\text{on}\quad\partial\Omega.$

This problem is a Sturm-Liouville eigenvalue problem ($-\Delta+V$ is a
diffusion-reaction operator, see, for example, [37, 14]) which has a countable
infinite set of positive eigenvalues ${E_{j}}\in{\mathbb{R}^{+}}$

$0<{E_{1}}<{E_{2}}\leq\cdots\leq{E_{j}}\leq\cdots$ (2.3)

with an associated set of orthonormal eigenfunctions ${u_{j}}$

$({u_{j}},{u_{k}})=\int_{\Omega}{{u_{j}}(x){u_{k}}}(x)\
\text{d}{\boldsymbol{x}}={\delta_{jk}},$ (2.4)

where $\delta_{jk}$ is the Kronecker delta which is equal to 1 when $j=k$ and
0 otherwise. The set of all the eigenvalues is the spectrum of the
Hamiltonian. We normalise the eigenfunctions in the $L^{2}$ space; hence, the
eigenfunctions are orthonormal under the scalar inner product. Let us define
two bilinear forms

$a(w,v)=\int_{\Omega}\nabla w\cdot\nabla v+Vwv\
\text{d}{\boldsymbol{x}}\quad\text{and}\quad b(w,v)=(w,v)=\int_{\Omega}wv\
\text{d}{\boldsymbol{x}},\quad\forall w,v\in H^{1}_{0}(\Omega),$ (2.5)

where $H^{1}_{0}(\Omega)$ is a Sobolev space with functions vanishing at the
boundary $\partial\Omega.$ Using this notation, the eigenfunctions are also
orthogonal with each other with respect to the energy’s inner product, i.e.,

$a(u_{j},u_{k})={E_{j}}({u_{j}},{u_{k}})={E_{j}}{\delta_{jk}}.$ (2.6)

We remark that these orthogonalities are critical in the development of the
proposed neural networks in Section 3.

###  2.2 Isogeometric Analysis

At the continuous level, the weak formulation for the eigenvalue problem (2.2)
is: Find all eigenvalues $E\in{\mathbb{R}^{+}}$ and eigenfunctions $u\in
H^{1}_{0}(\Omega)$ such that,

$a(w,u)=Eb(w,u),\quad\forall\ w\in H^{1}_{0}(\Omega),$ (2.7)

while at the discrete level, the isogeometric analysis (IGA, c.f., [22, 8, 12,
10]) for the eigenvalue problem (2.2) is: Find all eigenvalues
$E_{h}\in{\mathbb{R}^{+}}$ and eigenfunctions $u_{h}\in W_{h}$ such that,

$a(w_{h},u_{h})=E_{h}b(w_{h},u_{h}),\quad\forall\ w_{h}\in W_{h}(\Omega),$
(2.8)

where $W_{h}\subset H^{1}_{0}(\Omega)$ is the trial and test space spanned by
the B-spline or non-uniform rational basis spline (NURBS) basis functions [9].

In this paper, for the purpose of comparison with the state of the art, we
adopt the recently-developed soft isogeometric analysis (softIGA, c.f., [11,
29]). SoftIGA has a similar variational formulation as (2.8) which leads to
the matrix eigenvalue problem

$\mathbf{K}\mathbf{U}=E_{h}\mathbf{M}\mathbf{U},$ (2.9)

where
$\mathbf{K}_{jk}=a(\phi_{j},\phi_{k}),\mathbf{M}_{jk}=b(\phi_{j},\phi_{k}),$
and $\mathbf{U}$ contains the coefficients of the eigenfunction $u_{h}$, for
its representation in the linear combination of the B-spline basis functions.
For simplicity, the matrix $\mathbf{K}$ (although it contains a scaled mass)
is referred to as the stiffness matrix while the matrix $\mathbf{M}$ is
referred as the mass matrix, and $(E_{h},\mathbf{u}_{h})$ is the unknown
eigenpair. We refer to [11, 29] for details.

###  2.3 A Few Challenges on Numerical Computation

There are a few challenges in the numerical computation, such as using IGA or
softIGA, of the Schrödinger equation with random potentials. Firstly, there is
a non-uniqueness of solutions. The Schrödinger equation with random potentials
can have multiple solutions that correspond to the same (or very similar)
energy level. This non-uniqueness can make it difficult to accurately identify
the correct eigenstates, particularly when dealing with complex potential
landscapes. Secondly, it can be a multi-scale problem with large
sampling/discretisation errors. Random potentials can vary over multiple
length scales, which leads to inaccurate numerical simulations and can make it
challenging to choose an appropriate discretisation size and mesh. Thirdly,
the Schrödinger equation for many-body problems in multiple dimensions suffers
from the curse of dimensionality. The number of parameters needed to describe
the random potential can be very large. As the number of dimensions increases,
the computational cost of solving the Schrödinger equation can increase
exponentially.

Last but not least, the Schrödinger equation with random potentials is a
complex mathematical problem, and solving it numerically can be
computationally expensive. This is particularly true for large system sizes
and high disorder strengths, which require a large number of numerical
simulations. Moreover, solving the Schrödinger equation with a different
random potential requires constructing a new matrix eigenvalue problem (2.9),
which can be a time-consuming process. This can become particularly
challenging when solving the equation for a large number of random potentials.
In such cases, one potential solution is to train a neural network that takes
the random potential as input and produces the corresponding eigenstates as
outputs. For instance, to solve the Schrödinger equation with $N$ different
random potentials, one could use a classic method like softIGA and apply it
$N$ times, or alternatively, train a neural network using data from a subset
of cases and then use the well-trained neural network to solve the remaining
cases. This approach can save a significant amount of computational time,
particularly when $N$ is large. With these challenges in mind, we introduce
the following neural-network-based method as an alternative method to solve
the time-independent Schrödinger equations.

##  3 The Physics-Informed Neural Network

Deep neural networks are a type of machine learning algorithm inspired by the
structure and function of the human brain. They are composed of multiple
layers of artificial neurons, each performing a nonlinear transformation of
the input data. The output of each layer is fed as input to the next layer,
allowing the network to learn increasingly complex features and relationships
in the data.

Some common types of deep neural networks include convolutional neural
networks (CNNs) for image processing and computer vision tasks [28], long
short-term memory models (LSTMs) for sequential data processing [21], deep
belief networks (DBNs) for unsupervised learning [20], and generative
adversarial networks (GANs) for image and video synthesis [17]. Deep neural
networks have demonstrated impressive performance in many applications. In
this paper, we adopt the physics-informed neural networks (PINN) with a novel
loss function design to solve the Schrödinger equations.

###  3.1 PINN

The Physics-Informed Neural Networks (PINNs) are a type of feed-forward neural
networks (FNNs) that incorporate physics laws, typically in the form of
partial differential equations (PDEs), into their loss functions [30, 35].
This combination of neural networks and PDEs makes PINNs a powerful tool for
solving complex PDEs in various scientific and engineering applications. The
overall process of using PINNs to solve a PDE can be summarized as follows.
First, a neural network architecture is constructed. This involves defining
the structure of the network, including the number and activation functions of
its layers, as well as the number of neurons in each layer. Next, the loss
function is defined. The loss function quantifies the discrepancy between the
neural network’s predictions and the actual solution to the PDE. In the case
of PINNs, the loss function incorporates terms that enforce the PDE itself, as
well as any associated initial or boundary conditions. Once the neural network
and loss function are established, the network is trained using an
optimization algorithm, such as the Adam optimiser [24]. During the training
process, the neural network adjusts its weights and biases to minimise the
loss function, gradually improving its ability to approximate the PDE
solution. After training, the accuracy of the neural network’s predictions can
be evaluated by comparing them against calculations made by current numerical
methods, such as isogeometric analysis [8]. If the predictions are not
satisfactory, further refinement can be performed. This may involve modifying
the network architecture, adjusting the loss function, or selecting a
different optimization algorithm. The refined network is then retrained to
improve its performance. Our PINN implementation is summarised in Algorithm 1.

Typically, the loss function in PINN is defined as

$L=L_{\text{de}}+L_{\text{reg}},$ (3.10)

where $L_{\text{de}}$ specifies the loss associated with the PDE and
$L_{\text{reg}}$ specifies the loss associated with the regularisation such as
the boundary conditions. The loss function term $L_{\text{de}}$ involves
derivatives which are usually evaluated by the automatic differentiation by
Tensorflow [1] or by Pytorch [34]. In the following subsection, we present a
goal-oriented loss function for the best performance in solving the
Schrödinger equations with random potentials. A scale factor is assigned to
each regularisation term in $L_{\text{reg}}$. These scale factors are
considered hyperparameters and are tuned to emphasise certain constraints that
the network should adhere to in its prediction.

###  3.2 Loss Function Design

To demonstrate the main idea, we simplify our focus to the one-dimensional
case. We assume that the potential $V(x)$ is randomly distributed over $m$
regions over the domain $[0,1]$, and each region has a random value. Each set
of random values then characterises the individual differential eigenvalue
problem to be solved for Anderson localised states. We partition the interval
$[0,1]$ into $n+1$ equally spaced points $x_{j}=jh,j=0,1,\cdots,n$ where
$h=1/n$ is the grid/mesh size. We denote by
${\boldsymbol{u}}_{h}=(u_{h}(x_{0}),u_{h}(x_{1}),\cdots,u_{h}(x_{n}))^{T}$ as
a vector of the approximate values of the true solution $u(x)$ at each of
these coordinates $x_{j}$. Similarly, we denote by
${\boldsymbol{V}}=(V(x_{0}),V(x_{1}),\cdots,V(x_{n}))^{T}$ as a vector of the
approximate values of the potential $V(x)$. In our PINN model, we choose $n$
such that each region has exactly 5 nodes each for accuracy as well as for
training efficiency. To approximate the second derivative vector
$u^{\prime\prime}(x)$, we use the center finite difference method with an
accuracy order of six, i.e., the error is of $\mathcal{O}(h^{6})$.

We construct a PINN that splits into two separate networks: one that
calculates the eigenvalue $E$, and the other one produces the eigenfunction
$u(x)$ at certain nodes. These two networks are trained and optimised
simultaneously, with their outputs used in the same loss function described
below. We consider the combination of these two networks as one network and
apply it to discover the eigenstates. Figure 3.1 draws this model structure,
highlighting the split into separate networks. The network section finding the
eigenvalue has two hidden layers with $n$ nodes, while the network for finding
the eigenfunction has three hidden layers with $n$ nodes. This allows the
network to scale with the size of the mesh. Each layer in both networks uses
the ReLU activation function $f(x)=$ max($x,0$), and the loss minimisation is
done using the Adam optimiser with a learning rate of $5\times 10^{-4}$.

Figure 3.1: Neural network structure for the prediction of eigenstates of
Hamiltonians.

The driving force of this network towards a correct solution is the design of
its loss function, which is built of multiple terms relating to different
aspects of a correct solution. The neural network’s optimisation function
works to minimise the loss function, so we choose terms that are positive
errors for incorrect predictions of $E$ and $u(x)$, and are zero for correct
predictions. Thus, the minimisation of these terms leads the network to
produce correct predictions for $E$ and $u(x)$. During training, we choose to
input an initial vector of ones (of length $n$), which the network uses as a
basis to build its predictions of $E$ and $u(x)$ from.

The central term of our loss function is the term:

$L_{\text{de}}=\sqrt{\frac{1}{n+1}\sum_{i=0}^{n}\left(E_{h}u_{h}(x_{i})+u_{h}^{\prime\prime}(x_{j})-V(x_{j})u_{h}(x_{j})\right)^{2}}$
(3.11)

which is always non-negative and is zero for solutions $E$ and $u(x)$ that
satisfy equation 2.2.  To encourage the network to satisfy the boundary
conditions $u(0)=u(1)=0$, we create the term:

$L_{\text{bound}}=u_{h}^{2}(x_{0})+u_{h}^{2}(x_{n})$ (3.12)

which the minimisation of the loss function leads to
$u_{h}(x_{0})=u_{h}(x_{n})=0$. The major issue with using just these terms is
that the trivial solution $u(x)=0$ satisfies these terms (reducing them to
zero). Thus, we need the following term to penalise the approximate solution:

$L_{\text{norm}}=\left(\int_{0}^{1}u^{2}_{h}(x)\text{d}{\boldsymbol{x}}-1\right)^{2}.$
(3.13)

The minimisation of the loss function leads to approximate solutions of
$L^{2}$ norm 1. Integration is calculated using the mid-point rule and this
calculation is majorly local. In the PINN model developed in [23], this
condition was imposed using the terms like $1/u^{2}(x)$ and $1/E^{2}$. We
point out that since the eigenstates are localised, an eigenstate $u(x)=0$ for
most of the domain, leads to very large and even infinite values due to the
design of the term $1/u^{2}(x)$. Consequently, this leads to difficulties in
minimising the loss function and the overall neural network training.

These terms are sufficient to push the network to generate an admissible
eigenstate. However, the model often produces the same eigenstates, preventing
the discovery of higher states. To overcome this issue, once an eigenstate
$(E_{h},u_{h}(x))$ is discovered (using a patience condition described below),
we add the eigenstate $u_{h}(x)$ to a list of eigenstates $S$. Then, we stop
the training, recompile the network and begin training again, where we take
advantage of the known physical fact that the eigenfunctions of Hamiltonians
are orthogonal and incorporate the following term into the loss function:

$L_{\text{orth}}=\sum_{\hat{u}_{h}(x)\in
S}\left(\int_{0}^{1}\hat{u}_{h}(x)u_{h}(x)dx\right)^{2},$ (3.14)

where $u_{h}(x)$ is the network’s current state to be produced. We iteratively
repeat this process of adding eigenstates to $S$ and recompiling the network.
This iteration is automated and can be done an arbitrary number of times.
Sometimes the network will skip an eigenstate and converge to a higher one,
however, with enough iterations the network tends to find these states
eventually. To decide whether a prediction is admissible, we determine whether
$L_{\text{de}}$ is below a certain threshold, as well as a patience condition
as in [23] that checks if the absolute value of the change in $L_{\text{de}}$
and that of the eigenvalue is below another threshold. This allows the network
to converge to solutions with $L_{\text{de}}$ much lower than the chosen
threshold, given that it is converging fast enough.

The loss function we have built so far performs adequately on problems with
enough spread in the eigenvalues between different states. However, in
situations where the eigenstates have nearly identical eigenvalues, the
network is unable to distinguish between the states, converging to an
undesirable linear combination of the eigenfunctions. This occurs because
eigenspaces form vector spaces, so a linear combination of eigenfunctions of
the same eigenvalue $E$ will be another eigenfunction of eigenvalue $E$,
giving another eigenstate that would satisfy $L_{\text{de}}$ and all the
regularisation terms defined. When studying Hamiltonians with potentials
generated randomly, there is a disorder for Anderson localisation to occur. We
leverage this fact in our loss function with the following procedure.

We let the model run with the above loss function for a certain number of
epochs $q$, which is a hyperparameter. The goal is to set this parameter so
that at epoch $q$, the network has converged to a linear combination of the
eigenfunctions, for example in Figure 3.2.

Figure 3.2: Left: An example of a linear combination of eigenstates where
their eigenenergies are approximately equal. Right: The corresponding
potential $V(x)$ with $m=80$ uniform elements, where $V(x)$ was randomly
generated in each element according to the normal distribution with mean 1 and
standard deviation 0.3, scaled by $10^{6}$.

Since we know the eigenstates of the Hamiltonian are localised, each spike in
3.2 represents a separate localised state. We thus scan through the
eigenfunction and generate a list $K$ of the regions of the domain where the
norm of the eigenfunction is greater than some set hyperparameter, splitting
apart the regions where the function returns to zero to separate the localised
spikes. Then, in our iterative process of discovering eigenvalues, at each
step, we choose one of the intervals in $K$ (say, $[x_{a},x_{b}]$) and add the
following term to the loss function:

$L_{\text{loc}}=\left(\int_{0}^{1}u^{2}_{h}(x)\text{d}{\boldsymbol{x}}-\int_{x_{a}}^{x_{b}}u^{2}_{h}(x)\text{d}{\boldsymbol{x}}\right)$
(3.15)

which encourages the network to set the eigenfunction to zero outside of the
localisation interval $[x_{a},x_{b}]$. Equivalently, we could have added terms
like $L_{\text{end}}$ for each nodal point outside of $[x_{a},x_{b}]$. We
remark that this term is non-negative. With all these loss terms in mind, we
have

$L_{\text{reg}}=\alpha_{\text{bound}}L_{\text{bound}}+\alpha_{\text{norm}}L_{\text{norm}}+\alpha_{\text{orth}}L_{\text{orth}}+\alpha_{\text{loc}}L_{\text{loc}},$
(3.16)

which is added to $L_{\text{de}}$ to give the overall loss function $L$ in
(3.10). The $\alpha$ terms in $L_{\text{reg}}$ scale each term in the loss
function, providing emphasis to certain physical constraints. Through
hyperparameter tuning, we found that the following values provided adequate
predictions for the number of nodes we used (order of magnitudes $10^{2}$ to
$10^{3}$)

$\alpha_{\text{bound}}=n^{2},\quad\alpha_{\text{norm}}=n^{3},\quad\alpha_{\text{orth}}=n^{3},\quad\alpha_{\text{loc}}=n^{4}.$
(3.17)

The idea of having a large weight on the regularisation terms of the loss
function is to ensure the network prioritises minimising these terms before
working to a solution satisfying the differential equation (i.e. before
minimising $L_{\text{de}}$). This also penalises the network heavily for
straying its solution from the regularisations in $L_{\text{reg}}$, working to
project the parameter search space to a subspace adhering to the
regularisations. The method described above is presented in Algorithm 1. The
purpose of this algorithm is to discover the eigenstates for a given
potential. Throughout the paper, we use this loss function in our proposed
PINN with the structure shown in Figure 3.1.

Algorithm 1 Eigenstate Discovery Using PINN

1:Let $k$ be the epoch where $L_{\text{loc}}$ is added to the loss function

2:Let $j=1$

3:while $j\leq$ num_eigenstates do

4: Initiate the model

5: while training do

6: Calculate prediction of $E$ and $u(x)$

7: Compute $L_{\text{de}}$, $L_{\text{bound}}$, $L_{\text{norm}}$

8: Compute $L_{\text{orth}}$ using all stored eigenfunctions

9: if current epoch $=k$ then

10: Scan predicted eigenfunction for localisation sections

11: Choose the localisation section $S$ with greatest norm

12: Reset model weights

13: else if current epoch $>k$ then

14: Compute $L_{\text{loc}}$ across $S$

15: end if

16: Backpropagate and step

17: if patience condition and $L_{\text{de}}<$ threshold then

18: Store prediction of $E$ and $u(x)$ for future $L_{\text{orth}}$
calculation

19: $j=j+1$

20: break while

21: end if

22: end while

23:end while

##  4 Numerical Experiments

We test our PINN model on various Hamiltonians with randomly distributed
potential and demonstrate the network’s robustness against potential with
different distributions. In particular, we consider potentials with
distributions that cause the Hamiltonian to admit eigenstates with almost
identical eigenenergies. This is most prevalent in potentials distributed
according to the Binomial distribution. We test the network for different
numbers of randomly generated regions ($m$), and choose the number of nodes in
the mesh ($n$) such that each region contains 5 nodes (that is, we choose
$n=5m$).

The localised states usually occur in the region where the potential takes its
smallest value. When there are multiple disjoint regions at (or very close to)
this minimum value, each of these regions can admit a localised state, where
they all can have similar eigenenergy. Figure 4.3 shows the first four PINN
approximated eigenstates where the potential is shown in the right plot of
Figure 3.2. The network can distinguish between the first three eigenstates,
regardless of how similar their eigenvalues are.

We observe that if $u(x)$ satisfies equation 2.2, then $-u(x)$ also does.
$-u(x)$ also has the same $L^{2}$ norm as $u(x)$, so it satisfies all other
terms in our loss function. This is relevant as the network will choose to
converge to either $u(x)$ or $-u(x)$. For consistency, we plot all
eigenfunctions in their positive form.

Figure 4.3: The first four PINN approximated eigenstates with the potential
defined in the right plot of Figure 3.2. Here, ev refers to the eigenvalue
associated with the plotted eigenfunction. Figure 4.4: The potential $V(x)$
with $m=40$ uniform elements, where $V(x)$ was randomly generated in each
element according to the Bernoulli distribution with probability 0.5, scaled
by $10^{6}$ and shifted up by $10^{5}$.  Figure 4.5: The first four PINN
approximated eigenstates with the potential defined in Figure 4.4. Here, ev
refers to the eigenvalue associated with the plotted eigenfunction. Figure
4.6: The eigenvalue prediction with respect to the epoch (training process)
using the potential in Figure 4.4. The dotted lines represent the true
eigenvalues.

Figure 4.4 shows a potential with a Bernoulli distribution, while Figure 4.5
shows the first four eigenstates, the first three of which have almost
identical eigenvalues. We expect these eigenstates to localise in the troughs
of the potential, which is what we observe in the eigenfunction plots. The
states that localise in troughs of similar length tend to have almost the same
eigenvalue, which our model captures.

Figure 4.6 shows the model’s prediction of the eigenvalue as the training
process evolves. After 2000 epochs, the model adds the $L_{\text{loc}}$ term.
The model produces converging eigenvalues until around epoch 9000. At this
point, the model’s prediction has met the convergence criteria and saves the
eigenstate to use in $L_{\text{orth}}$ and it resets its weights to search for
the next eigenstate. In our implementation, we include a function that
generates a video to display the network’s eigenfunction prediction as the
training process evolves, which is available at our
GitHub111https://github.com/liamharcombe4/eigen-network page.

Table 1 displays the PINN eigenvalue errors where we used the reference
eigenvalues generated by the softIGA method with a very fine mesh (2000
nodes). We also test the model’s robustness against the number of nodes. We
note that the model achieves slightly better results for the potential
distribution according to the normal distribution compared to the uniform
distribution. However, the results for the Bernoulli distributed potentials
are much better, all with relative errors below 1%. This demonstrates the
success of our approach to the problem of eigenstates with similar eigenvalue.
Table 2 displays the eigenvalue prediction results for a single potential
generated by the Bernoulli distribution with a different number of mesh
elements. As the number of elements increases, there is a clear improvement in
the accuracy of the PINN-predicted eigenvalues. The convergence rate is of
order approximately 1 but the theoretical analysis of this convergence is
subject to future study.

Potential Distribution | Nodes ($n+1$) | Regions ($m$) | Eigenvalue Error ($\%$)
---|---|---|---
Normal | 400 | 80 | 0.863
Normal | 200 | 40 | 1.16
Normal | 100 | 20 | 0.123
Bernoulli | 200 | 40 | 0.351
Uniform | 200 | 40 | 1.92
Table 1: The relative eigenvalue errors of the PINN-approximated first eigenenergy of Hamiltonians with potentials generated according to various distributions. Nodes ($n+1$) | Eigenvalue Error ($\%$)
---|---
100 | 0.783
200 | 0.343
300 | 0.0512
Table 2: The relative eigenvalue errors of the PINN-approximated first
eigenenergy of Hamiltonian with a fixed Bernoulli distributed potential over
$m=20$ regions.

###  4.1 Loss Terms Analysis

Another potential based on the Bernoulli distribution is plotted in Figure
4.7, with the network’s eigenstate predictions in Figure 4.8. During the
training process of the network for these approximations, the value of each
term of the loss function is saved every 100 epochs. Figure 4.9 displays the
evolution of these terms throughout the training process.

Figure 4.7: The potential $V(x)$ with $m=40$ uniform elements, where $V(x)$
was randomly generated in each element according to the Bernoulli distribution
with probability 0.5, scaled by $10^{6}$ and shifted up by $10^{5}$. Figure
4.8: The first four PINN approximated eigenstates with the potential defined
in Figure 4.7. Figure 4.9: The values of the individual terms in the loss
function as the training process evolves.

The spikes in Figure 4.9 represent places in the training process where either
the $L_{\text{loc}}$ term is added to the loss function or the network’s
parameters are reset after successfully finding an eigenstate and preparing to
search for another. We highlight that the regularisation terms are minimised
almost instantly, shown by how sharp the spikes are in their plots over
training time. This demonstrates that our scaling of the regularisation terms
is successfully restricting the network to a search space of solutions that
satisfy these regularisation terms, before working to minimise
$L_{\text{de}}$.

###  4.2 Comparison with SoftIGA

Figure 4.11 plots the PINN’s eigenstate prediction compared to the SoftIGA
calculation for the case with a potential plotted in 4.10. We can see the
PINN’s eigenfunction matching well the ones calculated by SoftIGA, and the
eigenvalues each with percentage errors less than 0.5% from the SoftIGA
eigenvalue (which is the same for all four of these eigenstates). However,
Figure 4.12 shows a significant difference between the PINN and SoftIGA
calculation of the fifth and sixth eigenstates for this potential. The PINN’s
eigenvalue predictions are within 2% of the SoftIGA calculations, however, the
SoftIGA eigenfunctions are not localised to single regions. Instead, the
eigenfunctions appear to be linear combinations of localised states, while the
PINN’s calculations are individual localised states. We believe that the
individual localised states that make up the linear combination that SoftIGA
finds, and the localised states that the PINN finds are all within the same
eigenspace (have the same eigenvalue), and thus their linear combinations are
also eigenstates with this eigenvalue. We believe this eigenspace is the span
of all the localised states in the troughs of the potential that are one
segment wide (for example, the intervals $[0,0.02]$, $[0.08,0.1]$,
$[0.26,0.28]$, $[0.38,0.4]$, etc., in Figure 4.10). We also believe that the
eigenstates in Figure 4.11 are all part of the eigenspace spanned by the
localised states in the troughs of the potential that are two segments wide
($[0.16,0.2]$, $[0.48,0.52]$, $[0.78,0.82]$, $[0.94,0.98]$). Since our PINN
works by scanning through each segment to specifically detect a localised
state, the network does not pick up linear combinations of separate localised
states, such as the SoftIGA calculations in Figure 4.12. Thus, our network is
capable of separating these eigenspaces into their localised basis elements,
while SoftIGA is not, displaying an area where PINNs can outperform
traditional computational methods.

Figure 4.10: The potential $V(x)$ with $m=40$ uniform elements using $n=200$
nodes, where $V(x)$ was randomly generated in each element according to the
Bernoulli distribution with probability 0.5, scaled by $10^{6}$ and shifted up
by $10^{5}$. Figure 4.11: The first four PINN predicted eigenstates compared
to the SoftIGA calculated states with the potential defined in Figure 4.10.
Here, ev refers to the eigenvalue calculated by the PINN for each state, with
the SoftIGA calculated eigenvalue of 103679.96 for each of these four states.
Figure 4.12: The PINN and SoftIGA predicted sixth and seventh eigenstates for
the potential in 4.10 with their respective eigenvalue calculation (ev).

###  4.3 An Example in Two Dimensions

In 2D, we consider a special case where the potential can be decomposed as
$V(x,y)=V_{1}(x)+V_{2}(y)$ in a rectangular domain $\Omega=[0,1]^{2}$. The
Schrödinger operator can be rewritten as

$-\Delta+V(x,y)=(-\partial_{xx}+V_{1}(x))+(-\partial_{yy}+V_{2}(y)),$

which allows the decomposition in the weak solution, leading to an equivalence
of solving two 1D eigenvalue problems. In particular, in the finite and
isogeometric element setting with tensor-product meshes, the 2D matrix
eigenvalue problem can be also decomposed into tensor-product structures,
leading to 1D solvers. In this setting, a base function in 2D can be rewritten
as a product, i.e., $\phi(x,y)=\phi_{1}(x)\phi_{2}(y)$. Following the
derivations in [5, 3], for the eigenvalue problem (2.8) in 2D, let
$W_{h}=\text{span}\\{\phi_{1,j}(x)\phi_{2,k}(y)\\}$. Then, the bilinear forms
can be decomposed as

$\displaystyle a(\phi_{1,k}(x)\phi_{2,k}(y),\phi_{1,j}(x)\phi_{2,j}(y))$
$\displaystyle=a_{x}(\phi_{1,k}(x),\phi_{1,j}(x))\cdot
b_{y}(\phi_{2,k}(y),\phi_{2,j}(y))$ (4.18)
$\displaystyle\quad+b_{x}(\phi_{1,k}(x),\phi_{1,j}(x))\cdot
a_{y}(\phi_{2,k}(y),\phi_{2,j}(y)),$ $\displaystyle
b(\phi_{1,k}(x)\phi_{2,k}(y),\phi_{1,j}(x)\phi_{2,j}(y))$
$\displaystyle=b_{x}(\phi_{1,k}(x),\phi_{1,j}(x))\cdot
b_{y}(\phi_{2,k}(y),\phi_{2,j}(y)),$

where the 1D bilinear forms are defined as

$\displaystyle a_{x}(\phi_{1,k}(x),\phi_{1,j}(x))$
$\displaystyle=\int_{0}^{1}\Big{(}\frac{\text{d}}{\text{d}x}\phi_{1,k}(x)\frac{\text{d}}{\text{d}x}\phi_{1,j}(x)+V_{1}(x)\phi_{1,k}(x)\phi_{1,j}(x)\Big{)}\
\text{d}x,$ (4.19) $\displaystyle\quad a_{y}(\phi_{2,k}(y),\phi_{2,j}(y))$
$\displaystyle=\int_{0}^{1}\Big{(}\frac{\text{d}}{\text{d}y}\phi_{2,k}(y)\frac{\text{d}}{\text{d}y}\phi_{2,j}(y)+V_{2}(x)\phi_{2,k}(y)\phi_{2,j}(y)\Big{)}\
\text{d}y,$ $\displaystyle b_{x}(\phi_{1,k}(x),\phi_{1,j}(x))$
$\displaystyle=\int_{0}^{1}\phi_{1,k}(x)\phi_{1,j}(x)\ \text{d}x,$
$\displaystyle b_{y}(\phi_{2,k}(y),\phi_{2,j}(y))$
$\displaystyle=\int_{0}^{1}\phi_{2,k}(y)\phi_{2,j}(y)\ \text{d}y.$

This is also referred to as variational separability. With this in mind, we
can rewrite the matrix eigenvalue problem (2.9) in 2D as a problem with 1D
tensor-products:

$(K_{x}\otimes M_{y}+M_{x}\otimes K_{y})U=E_{h}(M_{x}\otimes M_{y})U,$ (4.20)

which allows

$\big{(}(K_{x}-E_{x}M_{x})\otimes
M_{y}+M_{x}\otimes(K_{y}-E_{y}M_{y})\big{)}U=0,$ (4.21)

where $E_{h}=E_{h,x}+E_{h,y}$. As the dimensions are separated, this leads to
two 1D eigenvalue problems (see [2, 5] for more details).

$\displaystyle K_{x}U_{x}$ $\displaystyle=E_{x}M_{x}U_{x},$ (4.22)
$\displaystyle K_{y}U_{y}$ $\displaystyle=E_{y}M_{y}U_{y}.$

With this insight of dimension decomposition in mind, we apply the PINN in $x$
and $y$ dimensions to find eigenstates for the 2D problem. The eigenvalues in
2D are the sum of the 1D eigenvalues while the eigenstates in 2D are the
tensor-products of the 1D eigenstates. This significantly reduces the
computational costs. As an example, Figure 4.13 shows the surface plot of one
of these two-dimensional potentials, with the PINN eigenstates displayed in
Figure 4.14. When the two 1D potentials admit localised eigenstates, the
eigenstates of the constructed 2D potential are also localised, as
$u(x,y)=u_{1}(x)u_{2}(y)$ is only nonzero on the small rectangular region
where both $u_{1}(x)$ and $u_{2}(y)$ are nonzero. This is demonstrated by
Figure 4.14, showing the eigenstates localised to the rectangular regions.

Figure 4.13: The surface plot of the potential $V(x,y)=V_{1}(x)+V_{2}(y)$
where $V_{1}$ and $V_{2}$ are separately randomly generated across $m=40$
regions using $n=200$ nodes according to the normal distribution with mean 1
and standard deviation 0.3, then scaled by $10^{6}$. Figure 4.14: The first
four PINN predicted eigenstates with the potential defined in Figure 4.13.

##  5 Concluding Remarks

Recently, there has been a growing interest in using machine learning and
neural networks to solve differential equations, particularly in the study of
Hamiltonians. In this paper, we extend current models by constructing a neural
network capable of distinguishing eigenstates from Hamiltonians with randomly
distributed potentials that lead to localisation in the states. Our approach
involves incorporating a normalisation loss term to avoid trivial solutions
and an orthogonality loss term to search for higher eigenstates. Our method
eliminates the need for driving terms, which can hinder a network’s ability to
converge to an admissible solution. The novelty of our approach lies in a new
localisation loss term, which encourages the network to converge to a solution
that is localised to a specific region. This idea is based on the physical
fact that localised eigenstates occur in randomly distributed potentials. We
utilise the network’s initial attempt at convergence to identify regions of
localised states before adding in this term to achieve the desired
approximation. We demonstrate the effectiveness of our approach in discovering
eigenstates for potentials generated randomly according to different
distributions, with the highlight being the Bernoulli distribution, which
leads to eigenstates with identical eigenvalue. Our model successfully
predicts these eigenstates with good accuracy, overcoming a challenge faced by
current models.

There are several potential avenues for future work. Some potential future
directions in this research area include: (1) exploring the use of more
complex loss functions that can capture additional physical constraints, such
as symmetries or conservation laws; (2) investigating the use of different
network architectures, such as convolutional neural networks or attention-
based models, to improve accuracy and scalability; (3) generalising the
approach to higher-dimensional systems or non-linear differential equations;
and (4) investigating the performance of the approach with larger and more
complex systems, including systems with a larger number of particles or more
complex interactions. Overall, the use of ANNs to solve differential equations
and approximate eigenstates of Hamiltonians has shown to be a promising
approach, and there are many potential directions for future research in this
area.

## References

  * [1] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al.  Tensorflow: a system for large-scale machine learning.  In Osdi, volume 16, pages 265–283. Savannah, GA, USA, 2016.
  * [2] M. Ainsworth and H. A. Wajid.  Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration.  SIAM Journal on Numerical Analysis, 48(1):346–371, 2010.
  * [3] P. Behnoudfar, V. M. Calo, Q. Deng, and P. D. Minev.  A variationally separable splitting for the generalized-$\alpha$ method for parabolic equations.  International Journal for Numerical Methods in Engineering, 121(5):828–841, 2020.
  * [4] S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis.  Physics-informed neural networks (PINNs) for fluid mechanics: A review.  Acta Mechanica Sinica, 37(12):1727–1738, 2021.
  * [5] V. Calo, Q. Deng, and V. Puzyrev.  Dispersion optimized quadratures for isogeometric analysis.  Journal of Computational and Applied Mathematics, 355:283–300, 2019\.
  * [6] Y. Chen, L. Lu, G. E. Karniadakis, and L. D. Negro.  Physics-informed neural networks for inverse problems in nano-optics and metamaterials.  Opt. Express, 28(8):11618–11633, Apr 2020.
  * [7] Y. Chen, M. Mancini, X. Zhu, and Z. Akata.  Semi-supervised and unsupervised deep visual learning: A survey, 2022\.
  * [8] J. A. Cottrell, T. J. Hughes, and Y. Bazilevs.  Isogeometric analysis: toward integration of CAD and FEA.  John Wiley & Sons, 2009.
  * [9] C. de Boor.  A practical guide to splines, volume 27 of Applied Mathematical Sciences.  Springer-Verlag, New York, revised edition, 2001.
  * [10] Q. Deng.  Isogeometric analysis of bound states of a quantum three-body problem in 1D.  In Computational Science – ICCS, pages 333–346, Cham, 2022. Springer International Publishing.
  * [11] Q. Deng, P. Behnoudfar, and V. M. Calo.  SoftIGA: Soft isogeometric analysis.  Computer Methods in Applied Mechanics and Engineering, 403:115705, 2023.
  * [12] Q. Deng, V. Puzyrev, and V. Calo.  Isogeometric spectral approximation for elliptic differential operators.  Journal of Computational Science, 2018.
  * [13] S. Desai, M. Mattheakis, H. Joy, P. Protopapas, and S. Roberts.  One-shot transfer learning of physics-informed neural networks, 2022.
  * [14] L. C. Evans.  Partial differential equations, volume 19 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, second edition, 2010.
  * [15] D. Finol, Y. Lu, V. Mahadevan, and A. Srivastava.  Deep convolutional neural networks for eigenvalue problems in mechanics.  International Journal for Numerical Methods in Engineering, 118, 01 2018.
  * [16] C. Flamant, P. Protopapas, and D. Sondak.  Solving differential equations using neural network solution bundles, 2020\.
  * [17] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial networks, 2014.
  * [18] L. Grubišić, M. Hajba, and D. Lacmanović.  Deep neural network model for approximating eigenmodes localized by a confining potential.  Entropy, 23(1), 2021.
  * [19] J. Han, A. Jentzen, and W. E.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510, aug 2018.
  * [20] G. Hinton, S. Osindero, and Y.-W. Teh.  A fast learning algorithm for deep belief nets.  Neural computation, 18:1527–54, 08 2006.
  * [21] S. Hochreiter and J. Schmidhuber.  Long short-term memory.  Neural computation, 9:1735–80, 12 1997.
  * [22] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  Comput. Methods Appl. Mech. Engrg., 194(39-41):4135–4195, 2005\.
  * [23] H. Jin, M. Mattheakis, and P. Protopapas.  Physics-informed neural networks for quantum eigenvalue problems.  In IJCNN at IEEE World Congress on Computational Intelligence, 2022\.
  * [24] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  International Conference on Learning Representations, 12 2014.
  * [25] F. Kotthoff, F. Pollmann, and G. De Tomasi.  Distinguishing an Anderson insulator from a many-body localized phase through space-time snapshots with neural networks.  Phys. Rev. B, 104:224307, Dec 2021.
  * [26] I. Lagaris, A. Likas, and D. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
  * [27] I. E. Lagaris, A. Likas, and D. I. Fotiadis.  Artificial neural network methods in quantum mechanics.  Computer Physics Communications, 104(1-3):1–14, 1997.
  * [28] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [29] D. Li and Q. Deng.  Soft isogeometric analysis of the bound states of a quantum three-body problem in 1D.  Journal of Computational Science, page 102032, 2023.
  * [30] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis.  DeepXDE: A deep learning library for solving differential equations.  SIAM review, 63(1):208–228, 2021.
  * [31] M. Magill, F. Qureshi, and H. W. de Haan.  Neural networks trained to solve differential equations learn general representations, 2018.
  * [32] T. Mano and T. Ohtsuki.  Machine learning the dynamics of quantum kicked rotor.  Annals of Physics, 435:168500, 2021.  Special Issue on Localisation 2020.
  * [33] M. Mattheakis, D. Sondak, A. S. Dogra, and P. Protopapas.  Hamiltonian neural networks for solving equations of motion.  Physical Review E, 105(6), jun 2022.
  * [34] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.  Automatic differentiation in pytorch.  2017\.
  * [35] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational physics, 378:686–707, 2019.
  * [36] J. Sirignano and K. Spiliopoulos.  Dgm: A deep learning algorithm for solving partial differential equations.  Journal of Computational Physics, 375:1339–1364, 2018.
  * [37] G. Strang and G. J. Fix.  An analysis of the finite element method, volume 212.  Prentice-Hall Englewood Cliffs, NJ, 1973.
  * [38] L. Yang, X. Meng, and G. E. Karniadakis.  B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data.  Journal of Computational Physics, 425:109913, 2021.


11institutetext: Jadavpur University, Kolkata, West Bengal, India
11email: {sksahu.etce.rs<EMAIL_ADDRESS>
# MULTI-MODAL MULTI-CLASS PARKINSON DISEASE CLASSIFICATION USING CNN and
DECISION LEVEL FUSION

Sushanta Kumar Sahu 11 0000-0002-4384-8939
Ananda S. Chowdhury 11 0000-0002-5799-3467 Sushanta Kumar Sahu
0000-0002-4384-8939
Ananda S. Chowdhury  0000-0002-5799-3467

###### Abstract

Parkinson’s disease (PD) is the second most common neurodegenerative disorder,
as reported by the World Health Organization (WHO). In this paper, we propose
a direct three-Class PD classification using two different modalities, namely,
MRI and DTI. The three classes used for classification are PD, Scans Without
Evidence of Dopamine Deficit (SWEDD) and Healthy Control (HC). We use white
matter (WM) and gray matter (GM) from the MRI and fractional anisotropy (FA)
and mean diffusivity (MD) from the DTI to achieve our goal. We train four
separate CNNs on the above four types of data. At the decision level, the
outputs of the four CNN models are fused with an optimal weighted average
fusion technique. We achieve an accuracy of 95.53% for the direct three-class
classification of PD, HC and SWEDD on the publicly available PPMI database.
Extensive comparisons including a series of ablation studies clearly
demonstrate the effectiveness of our proposed solution.

###### Keywords:

Parkinson’s disease (PD) Direct three-Class classification Multi-modal
DataDecision level fusion.

##  1 Introduction

Parkinson’s disease is the second most common neurological disorder that
affects movement and can cause tremors, stiffness, and difficulty with
coordination [1]. Early diagnosis of PD is important for effective treatment,
as there is currently no cure for the disease. However, diagnosis can be
challenging due to the variability of symptoms and lack of definitive
biomarkers. According to the World Health Organization (WHO), PD affects
approximately 1% of people aged 60 years and older worldwide. However there
are approximately 10% of clinically diagnosed patients with early stage PD who
exhibit normal dopaminergic functional scans. This class, which signifies a
medical condition distinct from PD, is known as Scans Without Evidence of
Dopamine Deficit (SWEDD) [2, 3]. As a result of the evolution of this new
class, difficulty of diagnosing PD has increased manifold, leading to a three-
class classification problem of PD vs. SWEDD vs. HC with class overlaps [3].

MRI, SPECT and PET are commonly used imaging techniques for PD diagnosis.
However, PET and SPECT are not preferred by doctors due to invasiveness and
cost [4]. DTI is a newer technique that measures water molecule movement to
analyze white matter microstructure which gets affected in PD. In the
literature, quite a few works were reported on PD classification based on
machine learning (ML) and deep learning (DL) models applied to neuroimaging
data. Salat et al. [5] found correlations between gray or white matter changes
and age using Voxel-based Morphometry (VBM). Adeli et al. [1] used a recursive
feature elimination approach for two-class classification with 81.9 %
accuracy. Cigdem et al. [6] proposed a total intracranial volume method with
93.7% accuracy. Singh et al. [2] presented a ML framework for three two-class
classifications. Chakraborty et al. [7] presented an DL model with 95.29%
accuracy. A DL-based ensemble learning technique was reported by [8] with
97.8% accuracy.

Recent research indicates that combining features from more than one imaging
modality can improve the classification accuracy. For example, Li et al. [9]
showed that combining DTI and MRI features improves the classification
accuracy in Alzheimer’s disease. The authors of [10] used MRI and DTI, but
only considered MD data from DTI for PD classification. They used a stacked
sparse auto-encoder to achieve better classification accuracy. In light of the
above findings, we anticipate that MRI and DTI can be effectively combined to
better analyze PD. To increase decision accuracy, we mention here a few
decision level fusion techniques. Majority voting technique is the most common
techniques used in late fusion [11]. This strategy, however, may not be
appropriate for multi-class classification applications. Single classifiers
work well on most subjects, but error rates are enhanced for some difficult-
to-classify subjects due to overlap across many categories. Instead of using a
majority voting strategy, a modified scheme called modulated rank averaging is
employed in [12]. We feel the accuracy may be enhanced even further by fine-
tuning the weights computed in the modulated rank averaging approach.

As a summary, we can say that there is a clear dearth of direct three-class PD
classification strategies and that too with multi-modal data. In this paper,
we present a direct three-class PD classification using CNNs and decision
level fusion. We investigate full potential of multimodal data i.e., FA & MD
from DTI and WM & GM from MRI. We use four CNNs to analyze these four types of
data. Outputs from all these four models are finally fused using an Optimal
Weighted Average Fusion (OWAF) technique at the decision level. Since
neuroimaging datasets are small, data augmentation is adopted to ensure proper
training with the CNNS [12]. We now summarize our contributions as below:

  1. 1.

We address a direct three-class classification task (PD, HC and SWEDD) for
Parkinson’s disease, which is certainly more challenging than the current
trend of a single binary classification (for 2-class problem) or multiple
binary classifications (for 3-class problem).

  2. 2.

We make effective use of the underlying potential of multi-modal neuroimaging,
namely T1-weighted MRI and DTI. In particular, we train four different CNNs on
WM, GM data from MRI and FA, MD data from DTI. Such in-depth analysis of
multi-modal neuroimaging data is largely missing in the analysis of PD.

  3. 3.

Finally, at the decision level, the outputs of each CNN model are fused using
an Optimal Weighted Average Fusion (OWAF) strategy to achieve state-of-the-art
classification accuracy.

##  2 Proposed method

Our solution pipeline for an end-to-end direct three-class classification of
PD from DTI and MRI consists of four CNN networks. Each CNN network yields a
$3\times 1$ probability vector, which represents the probability that the data
falls into one of three classes, i.e., PD, HC or SWEDD. The probability
vectors are then combined using the OWAF technique. In section 2.1, we discuss
how WM and GM are obtained from MRI data and MD and FA are used from DTI data.
Section 2.2 describes ADASYN, an oversampling strategy. Section 2.3 presents
the proposed CNN architecture. In section 2.4, we discuss the decision level
fusion. Figure 1 illustrates the overall pipeline of our solution.

Figure 1: Direct three-class Parkinson’s disease classification framework

###  2.1 Data Pre-processing

In this work, voxel-based morphometry (VBM) is used to prepare MRI data. The
data is preprocessed using SPM-12 tools and images are normalized using the
diffeomorphic anatomical registration with exponentiated lie algebra (DARTEL)
method [13]. This SPM-12 tool segments the whole MRI data into GM, WM and
cerebrospinal fluid, as well as the anatomical normalization of all images to
the same stereotactic space employing linear affine translation, non-linear
warping, smoothing and statistical analysis. After registration, GM and WM
volumetric images were obtained and the unmodulated image is defined as the
density map of grey matter (GMD) and white matter (WMD). The PPMI database
contains all information regarding DTI indices, including FA and MD. Brain
scans from PD groups as well as HC have distinct voxel MD values. MD and FA
can be expressed mathematically expressed as [14]:

$MD=\frac{\lambda_{1}+\lambda_{2}+\lambda_{3}}{3}=\frac{D_{xx}+D_{yy}+D_{zz}}{3}$
(1)
$FA=\sqrt{\frac{1}{2}}\sqrt{\frac{{(\lambda_{1}-\lambda_{2})}^{2}+{(\lambda_{2}-\lambda_{3})}^{2}+{(\lambda_{3}-\lambda_{1})}^{2}}{{\lambda_{1}}^{2}+{\lambda_{2}}^{2}+{\lambda_{3}}^{2}}}$
(2)

In equations 1 the diagonal terms of the diffusion tensor are
$D_{xx},D_{yy}andD_{zz}$ and the sum of these diagonal terms constitutes its
trace. After prepossessing, four types of data are made available, namely,
grey matter (GM), white matter (WM), fractional anisotropy (FA) and mean
diffusivity (MD).

Figure 2: Dataset division strategy

###  2.2 Data balancing with ADASYN

We find each of GM, WM, FA, MD to be highly imbalanced across the three
classes. Further, the number of training samples required to feed a DL model
is insufficient. So, ADASYN, an oversampling method is used to increase the
number of samples for each minority class [15]. The primary idea behind the
use of an ADASYN technique is to compute the weighted distribution of minority
samples based on a wide range of out-of-elegance neighbors. A difficult-to-
train minority instance surrounded by more out-of-class examples is given a
better chance of being augmented through producing synthetic samples. Using a
set of pseudo-probabilistic rules, a predetermined number of instances are
generated for every minority class depending on the weighted distribution of
its neighbors. Following the implementation of this up-sampling approach, the
total number of samples in each of the classes is 141 volumetric images. The
details of data set division strategies are shown in Figure 2.

Let a dataset consists of $m$ samples of the form $({x_{i},y_{i}})$; where
$i\in[1,m]$, $x_{i}$ being the $i^{th}$ sample of the n-dimensional feature
space X and $y_{i}$ being the label (class) of the sample ${x_{i}}$. Let
$m_{min}$ and $m_{maj}$ be number of the minority and majority class samples
respectively. Then the ratio of minority to majority sample will be expressed
as: d = $m_{min}$ /$m_{maj}$. Let, $\beta$ be the balance level of the
synthetic samples. So, $\beta=1$ represents that both the classes are balanced
using ADASYN. Number of synthetic minority data to be created is given by: 3.

$G=(m_{maj}-m_{min})\times\beta$ (3)

For $\beta\neq 1$, synthetic samples will be created for each set of minority
data based on the Euclidean distance of their k-nearest neighbors. The
dominates of the majority class in each neighborhood is expressed as:
$r_{i}=\frac{\Delta_{i}}{k}$. where $\Delta_{i}$ is the number of examples in
the k-nearest neighbors of $x_{i}$ that belong to the majority class. Higher
value of ${r_{i}}$ in a neighborhood indicates more examples of the majority
class which makes them harder to learn. We next determine how many synthetic
samples per neighborhood need to be generated as: $G_{i}=G\times\hat{r_{i}}$.
Here $\hat{r_{i}}$ represent normalized version of ${r_{i}}$. This captures
the adaptive nature of the ADASYN algorithm, which means more data is created
for harder to learn neighbourhoods. We generate $G_{i}$ number of synthetic
data ${s_{i,j},j=1,2,\cdots,G_{i}}$ for each neighbourhood $i$ as shown in
equation 4 below:

${s_{i,j}}=x_{i}+(x_{zi}-x_{i})\times\lambda$ (4)

Where $x_{i}$ and $x_{zi}$ are two minority occurrences in the same
neighborhood and $\lambda$ is a random integer between 0 and 1. In our work,
these synthetic images are generated with the help of the Scikit-Learn library
such that all the classes are balanced in nature.

###  2.3 Proposed CNN Architecture

We use four CNN models, each with ten convolutional layers and four dense (FC)
layers, for the direct three-class classification task of PD. The proposed
network’s architecture is depicted in Figure 3. We chose fewer parameters in
the proposed architecture than in the original VGG16 [16] by decreasing
network depth. Our proposed network is similar to that proposed in [17]. The
number of layers are limited to maintain a trade-off between accuracy and
computational cost. To generate feature representations of brain MR scans,
convolution layers are used. The final FC layer and a soft max operation are
used for the classification task. The volumetric input data is processed
slice-wise, with each slice having a size of $176\times 176$ pixels. We have
used max-pooling in the pooling layers to reduce the image size. The flattened
layers convert the reduced feature maps to a one-dimensional feature map. The
fully connected layers classify this feature map into three classes: HC, PD
and SWEDD. The cross-entropy loss function (CELF) is the most common loss
function used for classification problems since it has better convergence
speeds for training deep CNNs than MSE and hinge loss. As a result, we
consider CELF for this work, which is mathematically expressed as:

$L_{CELF}=-\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{K}p(y_{j}|\mathbf{x})\log\widehat{p}\left({y_{j}|\mathbf{x}}\right)$
(5)

where $p(y_{j}|\mathbf{x})$ is the original class label distribution and
$\widehat{p}(y_{j}|\mathbf{x})$ is the predicted label distribution from the
CNN network. Here $N$ and $K$ represent the number of models and number of
classes respectively. For our problem, $N=4$ and $K=3$.

Figure 3: Architecture of the proposed CNN

###  2.4 Decision Level Fusion of CNN Networks

We use four CNN models, one each for GM, WM, FA and MD. We fuse these
predicted probabilities with the help of suitable weights. The weights are
generated in two stages. In the first stage, the weights are generated using
the modulated rank averaging (MRA) method [12]. The weights in the MRA method
are given by:

$w_{i}=\frac{f_{i}}{\sum_{i=1}^{N-1}f_{i}+R_{max}}$ (6)

In equation 6, $f_{i}$ and $R_{max}$ indicate the normalizing factor and the
rank of the model having highest accuracy respectively. The normalizing factor
is calculated based on the rank of the current model and the difference
between the accuracy of the current and next model. In the second stage, these
weights are optimised using the grid search method [18]. Let the final
optimized weights be denoted by $w_{i}^{{}^{\prime}}$. Note that this weight
vector is fixed for all 3 classes. We combine this optimal weight vector with
the respective probability vectors to obtain the overall probability of
occurrence of respective classes. Let us denote by $PF_{j};j=1,2,3$, the
overall probability of occurrence of the $j^{th}$ class as a result of fusion.

$PF_{j}=\sum_{i=1}^{4}w_{i}^{{}^{\prime}}\times p_{ij}$ (7)

The final class will be the one for which $PF_{j}$ is maximum.

Table 1: Demographics of the subjects | Healthy Control | Parkinson’s Disease | SWEDD
---|---|---|---
Number of subjects | 67 | 177 | 37
Sex (Female/Male) | 24/43 | 65/122 | 14/23
Age* | 60.12 ±10.71 | 61.24 ±9.47 | 59.97 ±10.71
*mean ± standard deviation

##  3 Experimental Results

In this section we present the experimental results for three class PD
classification. The section is divided into three subsections. In the first
subsection, we provide an overview of the PPMI database. We then discuss data
preparation, computing configuration, parameter settings of ADASYN and CNN and
evaluation metrics. The second subsection shows a series of ablation studies
to demonstrate separate impacts of both MRI and DTI data and the proposed
fusion strategy. Finally, we include comparisons with several state-of-the-art
methods to demonstrate the effectiveness of the proposed solution.

###  3.1 Data preparation and implementation details

In our study, we included 281 subjects with baseline visits having both DTI
and MRI data from PPMI. This includes 67 HC, 177 PD and 37 SWEDD subjects.
Table 1 shows the demographics of the individuals used in this investigation.
For all experiments, we use a system with 16.0 GB DDR4 RAM, an Intel® Core™
i7-10750H CPU @ 2.60GHz and GPU of NVIDIA GeForse RTX 3060 @ 6GB. Here, 80% of
the data were randomly selected from the PD, SWEDD and HC groups to produce
the training set from 225 volumetric image scans. The remaining 56 volumetric
image scans, representing approximately 20% of each class, were utilized to
create the test set. For ADASYN, we experiment with different neighbor counts
(k) on the training set while keeping other settings unaltered. We find that k
= 30 produced the best results. Figure 2 shows the details of data set
division strategies. As a result of this, the total number of volumetric
images becomes 423, where each of the three classes have 141 volumetric
images. For the four CNN models, one each on GM, WM, FA and MD data, we train
the network for 100 epochs. ADAM optimizer and ReLU activation functions are
employed. The learning rate is initialized at $1\times 10^{-4}$ with a batch
size of 32. We use the same training parameters for each model (WM, GM, FA and
MD) such that there are no conflicts when we combine the outputs of the models
and make a fusion at the decision level. We evaluate the classification
performance using four different metrics. These are accuracy, precision,
recall and F1 score. All the measures are calculated using the Scikit Learn
packages [19].

Table 2: Impact of DTI and MRI data on direct 3-class classification

Data | Accuracy in (%) | Precision in (%) | Recall in (%) | F1 Score in (%)
---|---|---|---|---
MRI (WM) | 88.6 | 86.57 | 82.50 | 88.94
MRI (GM) | 84.2 | 81.76 | 82.50 | 84.80
DTI (MD) | 88.2 | 88.85 | 89.50 | 88.11
DTI(FA) | 80.94 | 80.35 | 81.00 | 81.59
WM and GM (MRI) | 90.02 | 90.17 | 90.02 | 90.09
FA and MD (DTI) | 91.14 | 91.23 | 91.14 | 91.18
| WM, GM, FA, MD
---
(MRI + DTI) using OWAF
95.53 | 93.64 | 91.99 | 92.74

Table 3: Impact of fusion strategies on direct 3-class classification Fusion Technique |  | Accuracy in (%)
---
| Precision in (%)
---
| Recall in (%)
---
| F1 Measure in (%)
---
Majority Voting | 92.19 | 92.5 | 94.87 | 92.34
Model Average Fusion | 88.6 | 82.5 | 86.57 | 85.44
Modulated Rank Average | 94.6 | 91.5 | 93.92 | 93.02
| OWAF (proposed)
---
95.53 | 93.64 | 91.99 | 92.74

###  3.2 Ablation Studies

We include two ablation studies. The first study demonstrates the utility of
using both MRI and DTI data. The second study conveys the benefit of OWAF, the
proposed fusion strategy. The four CNNs are trained and evaluated on both
single and multi-modal data from MRI and DTI. Our goal is to investigate the
effects of using single and multi-modal data on direct three-class. The
results of direct three-class classification are presented in Table 2. This
tables clearly illustrate that use of both DTI and MRI yields superior results
as compared to using MRI and DTI in isolation. Note that the magnitude of
improvement from the use of multi-modal data is clearly more significant for
the more challenging three-class classification problem. The four CNNs are
combined using four different fusion strategies at the decision level, namely;
majority voting, model average fusion, modulated rank averaging [12] and the
proposed optimal weighted average fusion (OWAF) based on the grid search
approach. When voting techniques are applied, the universal decision rule is
established simply by fusing the decisions made by separate models. In the
model average fusion method, the output probabilities of each model are simply
multiplied by the weight provided to that model based on its accuracy. In the
modulated rank averaging method (MRA), the output probabilities of each model
are updated using a weight generated based on their rank and the difference in
probabilities between the models. This method gives better results than the
model average fusion method and the majority voting method. In this work, we
take the weights generated using MRA as the initial weights. These base
weights are further optimised using the grid search method. In grid search, we
fine tune the weights ($w_{i}$) in the range of $w_{i}\pm 0.05$ with a step
size of 0.01. The final weights are used in our OWAF technique. Table 3
illustrates the effects of various fusion strategies. For fare comparison, we
use both MRI and DTI data in all cases. The experimental results in the table
3 clearly reveal that the proposed OWAF outperforms other fusion strategies.

Table 4: Comparisons of the proposed method with State-of-the-art Approaches

Approach | ML/DL | MODALITY | PD vs HC | PD vs. SWEDD | HC vs. SWEDD | PD vs. HC vs SWEDD
---|---|---|---|---|---|---
Ac | Pr | Re | Ac | Ac | Ac
Adeli 2016 [1] | ML | MRI | 81.9 | - | - | - | - | -
Cigdem 2018 [6] | ML | MRI | 93.7 | - | 95 | - | - | -
Prashanth 2018 [20] | ML | SPECT | 95 | - | 96.7 | - | - | -
Singh 2018 [2] | ML | MRI | 95.37 | - | - | 96.04 | 93.03 | -
Gabriel 2021 [21] | ML | MRI |  | 99.01(M)
---
87.10(F)
93.05 (A)
| 100(M)
---
97.2(F)
-
| 99.3(M)
---
100(F)
-
- | - | -
Li 2019 [10] | DL (AE) | MRI + DTI | 85.24 | 95.8 | 68.1 | - | 89.67 | -
Tremblay 2020 [22] | DL | MRI | 88.3 | 88.2 | 88.4 | - | - | -
Chakraborty 2020 [7] | DL | MRI | 95.3 | 92.7 | 91.4 | - | - | -
Sivaranjini 2020 [23] | DL | MRI | 88.9 | - | 89.3 | - | - | -
Rajanbabu 2022 [8] | DL (EL) | MRI | 97.5 | 97.9 | 97.1 | - | - | -
Proposed method | DL | MRI + DTI | 97.8 | 97.2 | 97.6 | 94.5 | 95.7 |  | 95.53
---
Pr: 93.64
Re: 91.99

  * •

Ac, Pr, Re, M, F, A, AE, EL, - Indicates Accuracy, Precision, Recall, Male,
Female, Average, Auto Encoder, Ensemble Learning & data not available
respectively. All the values are in %.

###  3.3 Comparisons with State-of-the-art Approaches

We compare our method with ten state-of-the-art approaches. There are no
results available for a direct 3-class PD classification. So, we compare our
results with those papers that have addressed the PD classification on the
PPMI database using single or multiple modalities and with three or fewer two-
class classifications. The results of comparisons are shown in Table 4. Out of
the ten methods we have considered, five are based on machine learning (ML)
and the rest five are based on deep learning (DL). Further, in four out of
five DL based approaches, only a single modality, namely, MRI is used for
classification. Also note that eight of these ten techniques have only
addressed a single two-class classification problem between PD and HC and did
not consider the challenging SWEDD class at all. The remaining two approaches
did consider SWEDD as a third class but have divided the three-class
classification problem into multiple binary classes [2, 10]. However, Li at
al. [10] did not report the classification results for PD vs. SWEDD in their
paper. In order to have fair comparisons, we have also included three binary
classifications as obtained from our method in this table. Our direct three-
class classification accuracy turns out to be superior than two-class
classification accuracy of at least eight out of 10 methods. It is also higher
than two out of three binary classification accuracy of [2]. Note that in [2],
the authors used a somewhat different experimental protocol by considering two
publicly available databases of ADNI and PPMI. In our work, we explicitly
consider data with both MRI and DTI for the same individual as available
solely in the PPMI database. Though the authors in [21] reported superior
classification accuracy for male, the accuracy is much less for female and
also for the average case (both male and female taken into account). Only, [8]
has reported a higher classification accuracy than ours. But, they have
considered only a single binary class classification of PD vs HC and ignored
the more challenging SWEDD class. If we consider the binary classification
results of our method, then we straightway outperform nine of the ten state-
of-the-art competitors and even beat the remaining method [2] in two out of
three classifications.

##  4 Conclusion

In this paper, we presented an automated solution for the direct three-class
classification of PD using both MRI and DTI data. Four different CNNs were
used for separate classifications with WM and GM data from MRI and MD and FA
data from DTI. An optimal weighted average decision fusion method was applied
next to integrate the individual classification outcomes. An overall three-
class classification accuracy of 95.53% is achieved. Extensive testing,
including a number of ablation studies on the publicly available PPMI database
clearly establishes the efficacy of our proposed formulation.

In future, we plan to deploy our model in clinical practice with data from
other neuro-imaging modalities. We also plan to to extend our model for
classifying other Parkinsonian syndromes.

## References

  * [1] E. Adeli, F. Shi, L. An, C. Y. Wee, G. Wu, T. Wang, D. Shen, Joint feature-sample selection and robust diagnosis of parkinson’s disease from mri data, NeuroImage 141 (2016) 206–219.
  * [2] G. Singh, L. Samavedham, E. C.-H. Lim, Determination of imaging biomarkers to decipher disease trajectories and differential diagnosis of neurodegenerative diseases (disease trend), Journal of Neuroscience Methods 305 (2018) 105–116.
  * [3] M. Kim, H. Park, Using tractography to distinguish swedd from parkinson’s disease patients based on connectivity, Parkinson’s Disease 2016, article id: 8704910 (2016).
  * [4] D. Long, J. Wang, M. Xuan, Q. Gu, X. Xu, D. Kong, M. Zhang, Automatic classification of early parkinson’s disease with multi-modal mr imaging, PLOS ONE 7 (2012) e47714.
  * [5] D. H. Salat, S. Y. Lee, A. Van der Kouwe, D. N. Greve, B. Fischl, H. D. Rosas, Age-associated alterations in cortical gray and white matter signal intensity and gray to white matter contrast, Neuroimage 48 (1) (2009) 21–28.
  * [6] O. Cigdem, I. Beheshti, H. Demirel, Effects of different covariates and contrasts on classification of parkinson’s disease using structural mri, Computers in Biology and Medicine 99 (2018) 173–181.
  * [7] S. Chakraborty, S. Aich, H.-C. Kim, Detection of parkinson’s disease from 3t t1 weighted mri scans using 3d convolutional neural network, Diagnostics 10 (6) (2020) 402.
  * [8] K. Rajanbabu, I. K. Veetil, V. Sowmya, E. A. Gopalakrishnan, K. P. Soman, Ensemble of deep transfer learning models for parkinson’s disease classification, in: V. S. Reddy, V. K. Prasad, J. Wang, K. T. V. Reddy (Eds.), Soft Computing and Signal Processing, Springer Singapore, Singapore, 2022, pp. 135–143.
  * [9] M. Li, Y. Qin, F. Gao, W. Zhu, X. He, Discriminative analysis of multivariate features from structural mri and diffusion tensor images, Magnetic resonance imaging 32 (8) (2014) 1043–1051.
  * [10] S. Li, H. Lei, F. Zhou, J. Gardezi, B. Lei, Longitudinal and multi-modal data learning for parkinson’s disease diagnosis via stacked sparse auto-encoder, in: 2019 IEEE 16th International Symposium on Biomedical Imaging (ISBI 2019), IEEE, 2019, pp. 384–387.
  * [11] A. Daskalakis, D. Glotsos, S. Kostopoulos, D. Cavouras, G. Nikiforidis, A comparative study of individual and ensemble majority vote cdna microarray image segmentation schemes, originating from a spot-adjustable based restoration framework, computer methods and programs in biomedicine 95 (1) (2009) 72–88.
  * [12] A. De, A. S. Chowdhury, DTI based alzheimer’s disease classification with rank modulated fusion of cnns and random forest, Expert Systems with Applications 169 (2021) 114338.
  * [13] J. Ashburner, A fast diffeomorphic image registration algorithm, Neuroimage 38 (1) (2007) 95–113.
  * [14] A. L. Alexander, J. E. Lee, M. Lazar, A. S. Field, Diffusion tensor imaging of the brain, Neurotherapeutics 4 (3) (2007) 316–329.
  * [15] Y. Pristyanto, A. F. Nugraha, A. Dahlan, L. A. Wirasakti, A. Ahmad Zein, I. Pratama, Multiclass imbalanced handling using adasyn oversampling and stacking algorithm, in: 2022 16th International Conference on Ubiquitous Information Management and Communication (IMCOM), 2022, pp. 1–5.
  * [16] K. Simonyan, A. Zisserman, Very deep convolutional networks for large-scale image recognition, arXiv preprint arXiv:1409.1556 (2014).
  * [17] S.-H. Wang, Q. Zhou, M. Yang, Y.-D. Zhang, Advian: Alzheimer’s disease vgg-inspired attention network based on convolutional block attention module and multiple way data augmentation, Frontiers in Aging Neuroscience 13 (2021) 687456\.
  * [18] H. A. Fayed, A. F. Atiya, Speed up grid-search for parameter selection of support vector machines, Applied Soft Computing 80 (2019) 202–210.
  * [19] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: Machine learning in Python, Journal of Machine Learning Research 12 (2011) 2825–2830.
  * [20] R. Prashanth, S. D. Roy, Early detection of parkinson’s disease through patient questionnaire and predictive modelling, International journal of medical informatics 119 (2018) 75–87.
  * [21] G. Solana-Lavalle, R. Rosas-Romero, Classification of ppmi mri scans with voxel-based morphometry and machine learning to assist in the diagnosis of parkinson’s disease, Computer Methods and Programs in Biomedicine 198 (1 2021).
  * [22] C. Tremblay, J. Mei, J. Frasnelli, Olfactory bulb surroundings can help to distinguish parkinson’s disease from non-parkinsonian olfactory dysfunction, NeuroImage: Clinical 28 (2020) 102457.
  * [23] S. Sivaranjini, C. Sujatha, Deep learning based diagnosis of parkinson’s disease using convolutional neural network, Multimedia tools and applications 79 (21) (2020) 15467–15479.