###### Abstract

In this paper, our focus is on investigating the impact of cosmological
constant on relativistic quantum systems comprising spin-0 scalar particles.
Our analysis centers around the Klein-Gordon equation, and we obtain both
approximate and exact analytical solutions for spin-0 particles of the quantum
system. Afterwards, we explore quantum oscillator fields by considering the
Klein-Gordon oscillator within the same space-time characterized by a
cosmological constant. We obtain an approximate expression for the energy
eigenvalue of the oscillator fields. In fact, the energy spectrum in both
scenarios are examined and show the influences of the cosmological constant
and geometry’s topology. Our investigation is situated within the context of a
magnetic universe-a four-dimensional cosmological space-time recognized as the
Bonnor-Melvin universe.

Quantum dynamics of spin-0 particles in a cosmological space-time

Faizuddin<EMAIL_ADDRESS><EMAIL_ADDRESS>
Department of Physics, University of Science & Technology Meghalaya, Ri-Bhoi,
793101, India
Abdelmalek<EMAIL_ADDRESS>;
<EMAIL_ADDRESS>
Laboratory of theoretical and applied Physics, Echahid Cheikh Larbi Tebessi
University, Algeria

Keywords: Quantum fields in curved space-time; Relativistic wave equations,
Solutions of wave equations: bound-state, special functions

PACS: 03.65.Pm; 03.65.Ge; 02.30.Gp.

##  1 Introduction

Delving into the profound intricacies of the interwoven tapestry between
gravitational forces and the dynamics of quantum mechanical systems captures
the imagination. Albert Einstein’s revolutionary general theory of relativity
(GR) skillfully paints gravity as an intrinsic geometric feature of space-time
[1]. This conceptual framework unravels the mesmerizing correlation between
space-time curvature and the genesis of classical gravitational fields,
furnishing accurate forecasts for phenomena such as gravitational waves [2]
and black holes [3]. Concurrently, the robust paradigm of quantum mechanics
(QM) [4] casts light on the subtle behaviors of particles at the microscopic
scale.

The convergence of these two domains of physics opens a gateway to profound
insights into the fundamental nature of the universe. The triumph of quantum
field theory in deciphering subatomic particle interactions and unraveling the
origins of weak, strong, and electromagnetic forces [5] amplifies the
anticipation surrounding this interdisciplinary fusion. Yet, the enduring
quest for a unified theory—a theory of quantum gravity harmonizing general
relativity and quantum mechanics—faces persistent challenges and technical
intricacies, with recent strides [6, 7]. These challenges propel intense
scientific endeavors as researchers diligently work to bridge remaining gaps,
aspiring to unveil the foundational framework that unifies these two pivotal
pillars of modern physics.

Certainly, among the most enthralling theories in physics is the general
theory of relativity. Its portrayal of gravitational interaction through the
geometry of space-time has uncovered an array of unexpected phenomena,
including the Lense-Thirring effect [8, 9] exploring frame dragging, the
revelation of gravitational waves [2], and the captivating phenomenon of
gravitational lensing [10, 11].

In recent times, there has been a noteworthy surge in interest surrounding the
investigation of magnetic fields, propelled by the discovery of systems
boasting exceptionally potent fields, exemplified by magnetars [12, 13], and
occurrences observed in heavy ion collisions [14, 15, 16]. This surge prompts
a captivating question: How can the magnetic field seamlessly integrate into
the overarching framework of general relativity? This inquiry introduces a
thought-provoking dimension to the ongoing exploration of the universe’s
fundamental forces and their interconnected dynamics. Exploring the
integration of the magnetic field within the framework of general relativity
raises intriguing questions. Presently, diverse solutions to the Einstein-
Maxwell equations exist, including the Manko solution [17, 18], the Bonnor-
Melvin universe [19, 20], and a recently proposed solution [46] that
incorporates the cosmological constant into the Bonnor-Melvin framework.
Shifting our focus to the intersection of general relativity and quantum
physics, a critical consideration emerges: how these two theories may
interrelate or if such a connection is even relevant. Numerous works have
addressed this inquiry, predominantly relying on the Klein-Gordon and Dirac
equations within curved spacetimes [21, 22]. This exploration extends to
diverse scenarios, encompassing particles in Schwarzschild [23] and Kerr black
holes [24], cosmic string backgrounds [25, 26, 27], quantum oscillators [28,
29, 30, 31, 32, 33, 34, 35, 36], the Casimir effect [37, 38], and particles
within the Hartle-Thorne spacetime [39]. These investigations, alongside
others [40, 41, 42], have yielded compelling insights into how quantum systems
respond to the arbitrary geometries of spacetime. Consequently, an intriguing
avenue of study involves the examination of quantum particles within a
spacetime influenced by a magnetic field. For instance, in [43], Dirac
particles were explored in the Melvin metric, while our current work delves
into the study of spin-0 bosons within a magnetic universe incorporating a
cosmological constant, as proposed in [46].

The dynamics of various relativistic quantum systems hinge critically on the
Dirac oscillator (DO). As highlighted by Ito et al.[44] in previous
advancements in the dynamics of spin-1/2 particles with a linear trajectory,
they demonstrated that the non-relativistic limit of the system results in an
ordinary harmonic oscillator with a significant spin-orbit coupling term.
Notably, Moshinsky and Szczepaniak [45] determined that the mentioned DO could
be derived from the free Dirac equation by introducing an external linear
potential, achieved through a minimal replacement of the momentum operator
$\hat{p}\longrightarrow\hat{p}-im\omega\beta\hat{r}$. It is important to
emphasize that, in addition to the theoretical exploration of the DO, valuable
insights can be gained by delving into its physical interpretation—an aspect
undeniably crucial for comprehending many relevant applications.

The Melvin magnetic universe is an exact static solution of Einstein-Maxwell
equations in which Maxwell’s magnetic pressure is balanced by gravitational
attraction. In this cylindrical symmetric space-time, the magnetic field is
parallel to the regular axis of symmetry. The stability as well as other
essential properties of Melvin’s magnetic universe have been extensively
studied by a number of authors. The space-time metric in cylindrical
coordinates is given by

$ds^{2}=\Lambda^{2}(r)\Big{(}-dt^{2}+dr^{2}+dz^{2}\Big{)}+\frac{r^{2}}{\Lambda^{2}(r)}\,d\phi^{2},$
(1)

where $\Lambda(r)=\Big{(}1+K^{2}\,r^{2}\Big{)}$ with $K=\frac{B}{2}$ and $B$
is the magnetic field strength. This metric reflects the curvature of the
space-time, determined by the existence of a magnetic field in the z
direction. In the limit of zero magnetic field, $B\to 0$, this metric (1)
reduces to the Minkowski space

$ds^{2}=-dt^{2}+dr^{2}+dz^{2}+r^{2}\,d\phi^{2}\,.$ (2)

The Bonnor-Melvin universe space-time [46] is described by

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=-dt^{2}+dz^{2}+\frac{1}{2\,\Lambda}\,\Big{(}d\rho^{2}+\alpha^{2}\,\sin^{2}\rho\,d\phi^{2}\Big{)},$
(3)

where $\Lambda$ denotes the cosmological constant, and $\alpha$ represents a
constant of integration intricately connected to the magnetic field through
the following relation [46]:

$H(r)=\frac{\alpha}{\sqrt{2}}\,\sin\rho\,.$ (4)

The covariant ($g_{\mu\nu}$) and contravariant form ($g_{\mu\nu}$) of the
metric tensor for the space-time (3) are given by

$\displaystyle g_{\mu\nu}=\begin{pmatrix}-1&0&0&0\\\
0&\frac{1}{2\,\Lambda}&0&0\\\
0&0&\frac{\alpha^{2}\,\sin^{2}r}{2\,\Lambda}&0\\\ 0&0&0&1\end{pmatrix},\quad
g^{\mu\nu}=\begin{pmatrix}-1&0&0&0\\\ 0&2\,\Lambda&0&0\\\
0&0&\frac{2\,\Lambda}{\alpha^{2}\,\sin^{2}r}&0\\\ 0&0&0&1\end{pmatrix}.$ (5)

The determinant of the metric tensor for the space-time (4) is given by

$det\,(g_{\mu\nu})=g=-\frac{\alpha^{2}}{4\,\Lambda^{2}}\,\sin^{2}r\,.$ (6)

Our current comprehension of the physical world relies on two foundational
frameworks: general relativity, describing gravitational interactions, and
quantum field theory, encompassing all non-gravitational interactions, leading
to the widely accepted standard model of particle physics and its proposed
generalizations and extensions. Despite their remarkable success in predicting
phenomena within their respective domains, inherent disparities between these
theories suggest an inherent conflict, rendering them incompatible. This not
only poses logical challenges but also limits our predictive capabilities in
scenarios where both gravitational and quantum effects coalesce. Such
scenarios include the physics of black holes, early universe cosmology, and
the intricate geometry of space-time at the Planck scale. The quest for a new
theory of quantum gravity becomes imperative to reconcile the disparities
between quantum field theory and general relativity and address unresolved
issues. Despite over a century of extensive efforts, the ability to make
quantitative predictions for quantum gravitational effects remains confined to
effective field theory frameworks.

The investigation into gravitational interactions within a quantum mechanical
system has been an active area of research, as evidenced by numerous studies
[47, 48, 49, 50, 51]. The exploration of quantum mechanical phenomena in
curved space-time introduces a novel interplay between quantum matter and
gravitational forces in the micro-particle realm. Notably, recent research has
been inspired by the Dirac oscillator embedded in a cosmic string background,
leading to a wealth of studies such as the dynamics of the Dirac oscillator in
cosmic string space-time [52, 53, 49, 27, 54, 55], the Aharonov-Casher effect
on the Dirac oscillator [47, 56], and the examination of non-inertial effects
on the Dirac oscillator in the backdrop of cosmic string space-time [57, 58].

Our motivation is to investigate the relativistic quantum dynamics of scalar
particles within the context of a magnetic universe, specifically described by
the Bonnor-Melvin metric with a cosmological constant. We aim to solve the
Klein-Gordon wave equation analytically and obtain the exact and approximate
eigenvalue solutions for the scalar particles. Subsequently, we explore the
relativistic quantum oscillator, modeled by the Klein-Gordon oscillator,
within the same magnetic universe background. We derive the radial wave
equation and obtain the approximate eigenvalue solution for the oscillator
field and conduct a comprehensive analysis of the results. Our findings reveal
that the energy eigenvalues in both analyses are influenced by the topology of
the geometry, as characterized by the parameter $\alpha$, and the cosmological
constant $\Lambda$. Furthermore, in the quantum oscillator system, the energy
spectrum undergoes modifications due to the oscillator frequency $\omega$, in
addition to the aforementioned parameters. The structure of this paper is
summarized as: In Sect. 2, we investigate the relativistic quantum dynamics of
scalar particle via the Klein-Gordon equation in the background of magnetic
universe. In Sect. 2, we focus on the relativistic quantum oscillator in the
background of same magnetic universe. In Sect. 4, we present our results and
discussion. Throughout this paper, we use the system of units, where $c=1=G$.

##  2 Scalar Particles in Cosmological Universe Background

In this part, we study the relativistic wave equation via the Klein-Gordon
equation in the background of a rotating wormhole space-time. Therefore, the
relativistic wave equation describing the quantum system is given by [62]

$\displaystyle\Bigg{[}-\frac{1}{\sqrt{-g}}\,\partial_{\mu}\,\Big{(}\sqrt{-g}\,g^{\mu\nu}\,\partial_{\nu}\Big{)}+M^{2}\Bigg{]}\,\Psi=0,$
(7)

where $M$ is the rest mass of the particles, $g$ is the determinant of the
metric tensor $g_{\mu\nu}$ with its inverse $g^{\mu\nu}$.

Therefore, the KG-equation (7) in the space-time (4) becomes can be expressed
in the following differential equation:

$\displaystyle\Bigg{[}-\frac{d^{2}}{dt^{2}}+2\,\Lambda\,\Bigg{\\{}\frac{d^{2}}{dr^{2}}+\frac{1}{\tan
r}\,\frac{d}{dr}+\frac{1}{\alpha^{2}\,\sin^{2}r}\,\frac{d^{2}}{d\phi^{2}}\Bigg{\\}}+\frac{d^{2}}{dz^{2}}-M^{2}\Bigg{]}\,\Psi=0\,.$
(8)

In quantum mechanical system, the total wave function is always expressible in
terms of different variables. In our case, we choose the following ansatz of
the total wave function $\Psi$ in terms of different variable functions
$\psi(x)$ and $\Theta(\theta)$ as follows:

$\Psi=\exp(-i\,E\,t)\,\exp(i\,m\,\phi)\,\exp(i\,k\,z)\,\psi(r),$ (9)

where $E$ is the particles energy, $m=0,\pm\,1,\pm\,2,....$ are the
eigenvalues of the angular quantum number. We are able to express total wave
function in the above since the differential equation (8) is independent of
the time $t$, angular $\phi$ coordinates, and $z$ coordinate.

Substituting the total wave function (9) into the differential equation (14)
and after separating the variables, we obtain the differential equations for
$\psi(x)$ as follows:

$\psi^{\prime\prime}+\frac{1}{\tan
r}\,\psi^{\prime}+\Bigg{[}\frac{1}{2\,\Lambda}\Big{\\{}E^{2}-k^{2}-M^{2}\Big{\\}}-\frac{\iota^{2}}{\sin^{2}r}\Bigg{]}\,\psi=0,$
(10)

where $\iota=\frac{|m|}{|\alpha|}$.

###  2.1 Approximate Eigenvalue Solution

In this part, we solve the above equation taking an approximation of the
trigonometric function up to the first order, that is, $\sin r\approx r$,
$\tan r\approx r$. Therefore, the radial wave equation (10) reduces to the
following form:

$r^{2}\,\psi^{\prime\prime}+r\,\psi^{\prime}+\Big{(}\lambda^{2}\,r^{2}-\iota^{2}\Big{)}\,\psi=0,$
(11)

where $\lambda=\frac{1}{2\,\Lambda}\,(E^{2}-k^{2}-M^{2})$.

This differential equation is known as Bessel second-order equation [60, 61]
whose solution is well-known. Here, we are interested on the well-behaved
solution of the above Bessel equation at the origin. Therefore, the Bessel
function of first which shows regularity at the origin is given by

$\psi=c_{1}\,J_{\iota}(|\lambda|\,r),$ (12)

where $c_{1}$ is an arbitrary constant. The asymptotic form of the Bessel
function of the first kind is given by [60]

$J_{\iota}\propto\cos\Big{(}|\lambda|\,r-\frac{\pi}{2}\,\iota-\frac{\pi}{4}\Big{)}$
(13)

We aim to confine the motion of scalar particles within a region characterized
by a hard-wall confining potential. This confinement is particularly
significant as it provides an excellent approximation when investigating the
quantum properties of systems such as gas molecules and other particles that
are inherently constrained within a defined spatial domain. The hard-wall
confinement is defined by a condition specifying that at a certain axial
distance, $r=r_{0}$, the radial wave function $\psi$ becomes zero, i.e.,
$\psi(r=r_{0})=0$. This condition is commonly referred to as the Dirichlet
condition in the literature. The study of the hard-wall potential has proven
valuable in various contexts, including its examination in the presence of
rotational effects on the scalar field [27], studies involving the Klein-
Gordon oscillator subjected to the influence of linear topological defects
[63, 26], investigations into non-inertial effects on a non-relativistic Dirac
particle [64], examinations of a Dirac neutral particle analogous to a quantum
dot [65], studies on the harmonic oscillator within an elastic medium
featuring a spiral dislocation [66], and investigations into the behavior of
Dirac and Klein-Gordon oscillators in the presence of a global monopole [67].
This exploration of the hard-wall potential in diverse scenarios enriches our
understanding of its impact on quantum systems, providing insights into the
behavior of scalar particles subject to this form of confinement. Therefore,
at $r=r_{0}$, we have $\psi(r=r_{0})=0$ and using (13), we obtain the
following relation:

$\Big{(}|\lambda|\,r-\frac{\pi}{2}\,\iota-\frac{\pi}{4}\Big{)}=(2\,n+1)\,\frac{\pi}{2},$
(14)

where $n=0,1,2,3,...$.

Simplification of the above relation gives us the energy eigenvalue associated
with the mode $\\{n,m\\}$ given by

$E_{n,m}=\pm\,\sqrt{k^{2}+M^{2}+\frac{2\,\Lambda\,\pi^{2}}{r^{2}_{0}}\,\Big{(}n+\frac{|m|}{2\,|\alpha|}+\frac{3}{4}\Big{)}^{2}}\,.$
(15)

(a)

(b)

(c)

(d)

Figure 1: Energy spectrum $E_{m,n}$ (equation (15)). Here the parameters are
set as $k=M=1$ and $r_{0}=0.5$.

(e)

(f)

Figure 2: Energy spectrum $E_{n}$ (equation (19)). Here the parameters are set
as $k=M=1$.

Equation (15) represents the relativistic energy eigenvalue of spin-0 scalar
particles in the Bonnor-Melvin magnetic universe. Notably, the energy spectrum
is intricately shaped by the underlying geometry, as characterized by the
parameters $\alpha$ and $\Lambda$. Additionally, the energy spectrum exhibits
variations with respect to the angular quantum number $m$. A noteworthy
observation is that the energy levels of the scalar particles exhibit equal
spacing on both sides of $E=0$. This intriguing feature signifies that in the
system under consideration, particles and antiparticles possess equal
energies. This symmetry in energy levels adds a distinctive characteristic to
the behavior of scalar particles in the Bonnor-Melvin magnetic universe,
emphasizing the equilibrium between particles and their corresponding
antiparticles within this specific framework.

We have generated Figure 1 to illustrate the energy spectrum for various
values of the cosmological constant $\lambda$ (Fig. 1(a)), while keeping
$\alpha=0.5$ and $m=1$. Additionally, Figure 1(b) depicts the impact of the
topological parameter $\alpha$ on the energy spectrum, with $\Lambda=0.5$ and
$m=1$. In Figure 1(c), we explore the influence of the angular quantum number
$m$ on the energy spectrum, maintaining fixed values of $\alpha=0.5=\Lambda$.
Lastly, Figure 1(d) illustrates the variation in the energy spectrum with
changes in the radial quantum number $n$, with $\alpha=0.5$ and $m=1$. These
figures showcase the nature of the energy spectrum as these parameters
$(\Lambda,\alpha,m,n)$ increase, providing valuable insights into the behavior
of the system under consideration.

###  2.2 Exact Eigenvalue Solution

In this part, we make a transformation as $u=\cos r$ into the differential
equation (10) results the following form:

$(1-u^{2})\,\psi^{\prime\prime}(u)-2\,u\,\psi^{\prime}(u)+\Big{(}\lambda-\frac{\iota^{2}}{1-u^{2}}\Big{)}\,\psi(u)=0\,.$
(16)

Setting $\ell=\Big{(}\lambda+\frac{1}{4}\Big{)}^{\frac{1}{2}}-\frac{1}{2}$
which results $\ell\,(\ell+1)=\lambda$ into the above equation (20), we obtain

$(1-u^{2})\,\psi^{\prime\prime}(u)-2\,u\,\psi^{\prime}(u)+\Big{[}\ell\,(\ell+1)-\frac{\iota^{2}}{1-u^{2}}\Big{]}\,\psi(u)=0\,.$
(17)

The above differential equation is the well-know associated Legendre
polynomial [60, 61] whose solutions are called associated Legendre polynomials
given by

$\psi(u)=P^{\iota}_{\ell}(u)=(-1)^{\iota}\,(1-u^{2})^{\frac{\iota}{2}}\,\frac{d^{\iota}}{du^{\iota}}\,(P_{\ell}(u)),$
(18)

where $P_{\ell}(u)$ is the Legendre polynomial.

It is known that this associated Legendre polynomial can be expressed in terms
of hypergeometric function
${}_{2}F_{1}\Big{(}-\ell,\ell+1;1-\iota;\frac{1-u}{2}\Big{)}$ [60] and by
series expansion one can show that this hypergeometric function ${}_{2}F_{1}$
becomes a finite degree polynomial of degree $n$ provided $-\ell$ must be a
non-positive integer, that is $\ell=n$, where $n=0,1,2,3,...$. Thus, we can
write $\ell\,(\ell+1)=\lambda=n\,(n+1)$.

Simplification of the above relation $\lambda=n\,(n+1)$ results the following
expression of the energy eigenvalue of the system given by

$E_{n}=\pm\,\sqrt{M^{2}+k^{2}+2\,\Lambda\,n\,(n+1)}\,.$ (19)

Equation (19) represents the relativistic energy eigenvalue of spin-0 scalar
particles in the Bonnor-Melvin magnetic universe. We see that this energy
expression is only influenced by the cosmological constant $\Lambda$ and
exhibits variations with respect to radial quantum number $n$. Here also, we
see that the energy levels of the scalar particles are equally spacing on both
sides about $E=0$. This signifies that in the system under consideration,
particles and antiparticles possess equal energies.

We present Figure 2 to depict the energy spectrum for different values of the
cosmological constant $\lambda$ (Fig. 2(a)) and the radial quantum number $n$
(Fig. 2(b)). The figures illustrate the behavior of the energy spectrum as
these parameters $(\Lambda,n)$ increase, offering insights into the system’s
dynamics with varying cosmological constants and radial quantum numbers.

##  3 Quantum Oscillator Fields in Cosmological Universe Background

In this part, we study the quantum oscillator fields via the Klein-Gordon
oscillator in the background of a cosmological solution. The oscillator field
is studied by replacing the radial momentum operator $\frac{\partial}{\partial
r}\to\Big{(}\frac{\partial}{\partial r}+M\,\omega\,r\Big{)}$, where $\omega$
is the oscillator frequency [28, 29, 30, 31, 32, 33, 34, 35, 36]. Therefore,
the relativistic wave equation describing the oscillator field is given by

$\displaystyle\Big{[}\frac{1}{\sqrt{-g}}\,(\partial_{\mu}+M\,\omega\,X_{\mu})\,(\sqrt{-g}\,g^{\mu\nu})\,(\partial_{\nu}-M\,\omega\,X_{\nu})\Big{]}\,\Psi=M^{2}\,\Psi,$
(20)

where $M$ is the rest mass of the particles, $\omega$ is the oscillator
frequency, $X_{\mu}=(0,r,0,0)$ is the four-vector, $g$ is the determinant of
the metric tensor $g_{\mu\nu}$ with its inverse $g^{\mu\nu}$.

Explicitly writing the wave equation (20) in the space-time background (4), we
obtain the following differential equation:

$\displaystyle\Bigg{[}-\frac{d^{2}}{dt^{2}}+2\,\Lambda\,\Bigg{\\{}\frac{d^{2}}{dr^{2}}+\frac{1}{\tan
r}\,\frac{d}{dr}-M\,\omega-\frac{M\,\omega\,r}{\tan
r}-M^{2}\,\omega^{2}\,r^{2}+\frac{1}{\alpha^{2}\,\sin^{2}r}\,\frac{d^{2}}{d\phi^{2}}\Bigg{\\}}$
$\displaystyle+\frac{d^{2}}{dz^{2}}-M^{2}\Bigg{]}\,\Psi=0\,.$ (21)

Substituting the wave function ansatz (9) into the above differential equation
(21) results the following second-order differential equation form:

$\psi^{\prime\prime}+\frac{1}{\tan
r}\,\psi^{\prime}+\Bigg{[}\frac{(E^{2}-M^{2})}{2\,\Lambda}-M\,\omega-\frac{M\,\omega\,r}{\tan
r}-M^{2}\,\omega^{2}\,r^{2}-\frac{k^{2}}{2\,\Lambda}-\frac{\iota^{2}}{\sin^{2}r}\Bigg{]}\,\psi=0.$
(22)

(a)

(b)

(c)

(d)

Figure 3: The energy spectrum $E_{n,m}$ for the relation given in Equation 26,
where the parameters are set as $k=M=1$ and $\omega=0.5$.

We solve the above equation by taking an approximation of the trigonometric
function up to the first order, that is, $\sin r\approx r$, $\tan r\approx r$.
Therefore, the radial wave equation (22) reduces to the following form:

$\psi^{\prime\prime}+\frac{1}{r}\,\psi^{\prime}+\Bigg{(}\lambda^{2}-M^{2}\,\omega^{2}\,r^{2}-\frac{\iota^{2}}{r^{2}}\Bigg{)}\,\psi=0,$
(23)

where $\lambda^{2}=\frac{E^{2}-M^{2}-k^{2}}{2\,\Lambda}-2\,M\,\omega$.

Transforming to a new variable via $s=M\,\omega\,r^{2}$ into the above
equation (23) results the following equation form [68]:

$\psi^{\prime\prime}(s)+\frac{(c_{1}-\,c_{2}\,s)}{s\,(1-c_{3}\,s)}\,\psi^{\prime}(s)+\frac{\Big{(}-\xi_{1}\,s^{2}+\xi_{2}\,s-\xi_{3}\Big{)}}{s^{2}\,(1-c_{3}\,s)^{2}}\,\psi(s)=0,$
(24)

where $c_{1}=1$, $c_{2}=0=c_{3}$ and

$\xi_{1}=\frac{1}{4},\quad\xi_{2}=\frac{\lambda}{4\,M\,\omega},\quad\xi_{3}=\frac{\iota^{2}}{4}.$
(25)

Therefore, the energy eigenvalue expression is given by

$E_{n,m}=\sqrt{M^{2}+k^{2}+8\,M\,\omega\,\Lambda\,\Big{(}n+1+\frac{|m|}{2|\alpha|}\Big{)}},$
(26)

where $n=0,1,2,....$.

The corresponding radial wave function will be

$\psi(s)=\mathcal{N}\,s^{\frac{|m|}{2|\alpha|}}\,e^{-s/2}\,L^{(\frac{|m|}{|\alpha|})}_{n}(s),$
(27)

where $\mathcal{N}$ is the normalization constant.

Equation (26) characterizes the relativistic energy eigenvalue of oscillator
fields within the Bonnor-Melvin magnetic universe. Notably, the expression for
energy is influenced by the underlying geometry, as described by the parameter
$\alpha$, and the cosmological constant $\Lambda$. In addition, the energy
spectrum of the oscillator fields undergoes modifications introduced by the
oscillator frequency $\omega$ and reveals variations in relation to the
quantum numbers $(n,m)$. Here also, we see that the energy levels are equally
spacing on both sides about $E=0$.

We present Figure 3, illustrating the energy spectrum for different values of
the cosmological constant $\lambda$ (Fig. 3(a)), where $\alpha=0.5$ and $m=1$.
Additionally, Figure 3(b) demonstrates the impact of the topological parameter
$\alpha$ on the energy spectrum, with $\Lambda=0.5$ and $m=1$. Figure 3(c)
portrays the influence of the angular quantum number $m$ on the energy
spectrum while maintaining fixed values of $\alpha=0.5$ and $\Lambda=0.5$.
Lastly, Figure 3(d) illustrates the variation in the energy spectrum with
changes in the radial quantum number $n$, while $\alpha=0.5$ and $m=1$. The
figures reveal that the linear nature of the energy spectrum undergoes changes
with increasing values of these parameters $(\Lambda,\alpha,m,n)$, providing
valuable insights into the evolving dynamics of the oscillator fields under
consideration.

##  4 Conclusions

Several studies have explored exact solutions to the Einstein-Maxwell field
equations, both with and without a cosmological constant, across various
dimensions. These solutions typically arise from the electromagnetic field
generated along axial and azimuthal directions. The pioneering work of Melvin
introduced the first solution with a purely magnetic field, subsequently
generalized as the Bonnor-Melvin magnetic universe. Recently, Zofka [46]
extended this framework by incorporating a non-zero cosmological constant into
the Bonnor-Melvin magnetic universe which is of our interest in the context of
quantum mechanical systems. Investigating quantum mechanical phenomena within
this magnetic universe provides valuable insights into the dynamics of scalar
particles under the influences cosmological constant.

This study aimed to investigate the behavior of spin-0 scalar particles within
the background of the Bonnor-Melvin magnetic universe with a non-zero
cosmological constant. The focus is on understanding the effects on the total
wave function, which provides comprehensive information about the quantum
system under investigation. Through this exploration, we aimed to shed light
on the intricate dynamics of scalar particles in this unique magnetic universe
and contribute to our understanding of quantum mechanical phenomena in such
space-time configurations.

In Section 2, we derived the radial equation of the Klein-Gordon equation in
this Bonnor-Melvin magnetic universe with a cosmological constant. We then
solved this radial equation and presented the approximate energy eigenvalue
given in Equation (15) and the exact analytical solution given in Equation
(19). We observed that the energy spectrum in Equation (15) is influenced by
the topology of the geometry $\alpha$ and the cosmological constant
$\lambda>0$. Furthermore, this energy spectrum is altered by the angular
quantum number $m$. However, the energy spectrum given in Equation (10) is
independent of both the topological parameter $\alpha$ and the angular quantum
number $m$. Another interesting observation is that the energy levels are
equally spaced on either side of $E=0$, indicating that particles and
antiparticles possess equal energies for fixed quantum states $\\{n,m\\}$.

In Section 3, we further investigated relativistic quantum oscillator fields
described by the Klein-Gordon oscillator. We derived the radial equation and
obtained the approximate energy eigenvalue of the oscillator fields given in
Equation (26). Here, too, we observed that the energy spectrum is influenced
by the topology of the geometry $\alpha$, the cosmological constant
$\lambda>0$, and the oscillator frequency $\omega$. Furthermore, this energy
spectrum is altered by the angular quantum number $m$. Also the energy levels
are equally spaced on either side of $E=0$, indicating that particles and
antiparticles possess equal energies for fixed quantum states $\\{n,m\\}$.

Thus, we presented an interesting topic, quantum field in curved space,
specifically focused on a magnetic universe with a non-zero cosmological
constant. We observed that the results get modified by various parameters
involved in the space-time geometry. The connection between these two theories
(quantum mechanics and gravity) is still not yet known completely.

## Conflict of Interest

There is no conflict of interests.

## Funding Statement

No funding agency is associated with this manuscript.

## Data Availability Statement

No data are generated or analysed during this study.

## References

  * [1] A. Einstein, Annalen Phys. 49, 769 (1916).
  * [2] B. P. Abbott et al, Phys Rev Lett 116, 061102 (2016).
  * [3] K. Akiyama et al, Astrophys J Lett 875, L1 (2019).
  * [4] R. P. Feynman, and A. R. Hibbs, Quantum mechanics and path integrals, Courier Corporation(1965).
  * [5] M. D. Schwartz, Quantum field theory and the standard model, Cambridge University Press (2013).
  * [6] A. Ashtekar, and J. J. Stachel, Conceptual problems of quantum gravity, Birkhäuser (1991).
  * [7] L. Smolin, The trouble with physics: The rise of string theory, the fall of a science, and what comes next, Mariner Books (2007).
  * [8] L. I. Schiff, Phys. Rev. Lett.4, 215(1960).
  * [9] J. Lense and H. Thirring, Physikalische Zeitschrift, 19, 156(1918).
  * [10] A. Einstein, Science, 84, 506 (1936).
  * [11] S. Refsdal and J. Surdej, Reports on Progress in Physics, 57, 117 (1994).
  * [12] C. Thompson and R. C. Duncan, Monthly Notices of the Royal Astronomical Society, 275, 255 (1995).
  * [13] C. Kouveliotou et al., Nature, 393, 235 (1998).
  * [14] U. Gürsoy, D. Kharzeev, and K. Rajagopal, Phys. Rev. C, 89, 054905 (2014).
  * [15] A. Bzdak and V. Skokov, Phys. Lett. B 710, 171 (2012).
  * [16] V. Voronyuk, V. D. Toneev, W. Cassing, E. L. Bratkovskaya, V. P. Konchakovski, and S. A. Voloshin, Phys. Rev. C 83, 054911 (2011).
  * [17] T. Gutsunaev and V. Manko, Phys. Lett. A 123, 215 (1987).
  * [18] T. Gutsunaev and V. Manko, Phys. Lett. A 132, 85 (1988).
  * [19] W. B. Bonnor, Proc. Phys. Soc. Section A 67, 225 (1954).
  * [20] M. Melvin, Phys. Lett. 8, 65 (1964).
  * [21] L. Parker, Phys. Rev. Lett. 44, 1559 (1980).
  * [22] N. D. Birrell, N. D. Birrell, and P. Davies, Quantum fields in curved space, Cambridge University Press (1984).
  * [23] E. Elizalde, Phys. Rev. D 36, 1269 (1987).
  * [24] S. Chandrasekhar, Proc. Roy. Soc. London. A. Math. Phys. Sci. 349, 571 (1976).
  * [25] L. C. N. Santos and C. C. Barros Jr., Eur. Phys. J. C 77, 186 (2017).
  * [26] L. C. N. Santos and C. C. Barros Jr., Eur. Phys. J. C 78, 13 (2018).
  * [27] R. L. L. Vitória and K. Bakke, Eur. Phys. J. C 78, 175 (2018).
  * [28] F. Ahmed, Int. J. Mod. Phys. A 37, 2250186 (2022).
  * [29] F. Ahmed, Commun. Theor. Phys. 75, 025202 (2023).
  * [30] L. C. N. Santos, C. E. Mota, and C. C. Barros Jr., Adv. High Energy Phys. 2019, 2729352 (2019).
  * [31] Y. Yang, Z.-W. Long, Q.-K. Ran, H. Chen, Z.-L. Zhao, and C.-Y. Long, Int. J. Mod. Phys. A 36, 2150023 (2021).
  * [32] A. R. Soares, R. L. L. Vitória, and H. Aounallah, Eur. Phys. J. Plus 136, 966 (2021).
  * [33] A. Bouzenada, A. Boumali, Ann. Phys. 452, 169302 (2023).
  * [34] A. Bouzenada, A. Boumali, R. L. L. Vitoria, F. Ahmed and M. Al-Raeei, Nucl. Phys. B 994, 116288 (2023).
  * [35] A. Bouzenada, A. Boumali and F. Serdouk, Theor. Math. Phys,216, 1055 (2023).
  * [36] A. Bouzenada, A. Boumali, and E. O. Silva, Ann. Phys. 458, 169479 (2023).
  * [37] V. B. Bezerra, M. S. Cunha, L. F. F. Freitas, C. R. Muniz, and M. O. Tahim, Mod. Phys. Lett. A 32, 1750005 (2016).
  * [38] L. C. N. Santos and C. C. Barros Jr., Int. J. Mod. Phys. A 33, 1850122 (2018).
  * [39] E. O. Pinho and C. C. Barros Jr., Eur. Phys. J. C. 83, 745 (2023).
  * [40] P. Sedaghatnia, H. Hassanabadi, and F. Ahmed, Eur. Phys. J. C. 79, 541 (2019).
  * [41] A. Guvendi, S. Zare and H. Hassanabadi, Phys. Dark Univ. 38, 101133 (2022).
  * [42] R. L. L. Vitória and K. Bakke, Eur. Phys. J. Plus 133, 490 (2018).
  * [43] L. C. N. Santos and C. C. Barros Jr., Eur. Phys. J. C 76, 560 (2016).
  * [44] D. Itô, K. Mori, and E. W. Carriere, Il Nuovo Cimento A (1965-1970) 51, 1119 (1967).
  * [45] M. Moshinsky and A. P. Szczepaniak, J. Phys. A 22 (1989).
  * [46] M. Žofka, Phys. Rev. D 99, 044058 (2019).
  * [47] K. Bakke and C. Furtado, Ann. Phys. (N. Y.) 336, 489 (2013).
  * [48] K. Bakke, Gen. Relativ. Grav. 45, 1847 (2013).
  * [49] K. Bakke, Eur. Phys. J. Plus 127, 82 (2012).
  * [50] E. R. Figueiredo Medeiros and E. R. Bezerra de Mello, Eur. Phys. J. C 72, 2051 (2012).
  * [51] L. B. Castro, Eur. Phys. J. C 75, 287 (2015).
  * [52] J. Carvalho, C. Furtado and F. Moraes, Phys. Rev. A 84, 032109 (2011).
  * [53] H. Hassanabadi, A. Afshardoost and S. Zarrinkamar, Ann. Phys. (N. Y.) 356, 346 (2015).
  * [54] K. Bakke and C. Furtado, Ann. Phys. (N. Y.) 336, 489 (2013).
  * [55] K. Bakke, H. Mota, Eur. Phys. J. Plus 133, 409 (2018).
  * [56] H. F. Mota and K. Bakke, Phys. Rev. D 89, 027702 (2014).
  * [57] K. Bakke, Ann. Phys. (N. Y.) 346, 51 (2014).
  * [58] P. Strange and L. H. Ryder, Phys. Lett. A 380, 3465 (2016).
  * [59] A. Ronveaux, Heun’s Differential Equations, Oxford University Press, Oxford (1995).
  * [60] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, New York: Dover (1972).
  * [61] G. B. Arfken, H. J. Weber and F. E. Harris, Mathematical Methods for Physicists, Elsevier (2012).
  * [62] W. Greiner, Relativistic Quantum Mechanics. Wave Equations, Springer-Verlag, Berlin, Gemnay (2000).
  * [63] R. L. L. Vitoria, K. Bakke, Int. J. Mod. Phys. D 27, 1850005 (2018).
  * [64] K. Bakke, Ann. Phys. 346, 51 (2014).
  * [65] K. Bakke, Eur. Phys. J. B 85, 354 (2012).
  * [66] A. V. D. M. Maia, K. Bakke, Physica B 531, 213 (2018).
  * [67] E. A. F. Bragança, R. L. L. Vitória, H. Belich, E. R. Bezerra de Mello, Eur. Phys. J. C 80, 206 (2020).
  * [68] A. F. Nikiforov and V. B. Uvarov, Special Functions of Mathematical Physics, Birkhauser (1988).


# Locally orderless tensor networks for classifying two- and three-dimensional
medical images

Raghavendra Selvan 1,2<EMAIL_ADDRESS>Silas Ørting 1<EMAIL_ADDRESS>Erik B
Dam 1<EMAIL_ADDRESS>
1 Department of Computer Science, University of Copenhagen, Denmark
2 Department of Neuroscience, University of Copenhagen, Denmark

###### Abstract

Tensor networks are factorisations of high rank tensors into networks of lower
rank tensors and have primarily been used to analyse quantum many-body
problems. Tensor networks have seen a recent surge of interest in relation to
supervised learning tasks with a focus on image classification. In this work,
we improve upon the matrix product state (MPS) tensor networks that can
operate on one-dimensional vectors to be useful for working with 2D and 3D
medical images. We treat small image regions as orderless, squeeze their
spatial information into feature dimensions and then perform MPS operations on
these locally orderless regions. These local representations are then
aggregated in a hierarchical manner to retain global structure. The proposed
locally orderless tensor network (LoTeNet111Source code is available here:
https://github.com/raghavian/LoTeNet_pytorch/) is compared with relevant
methods on three datasets. The architecture of LoTeNet is fixed in all
experiments and we show it requires lesser computational resources to attain
performance on par or superior to the compared methods.

Keywords: Tensor networks, Image classification, histopathology, thoracic CT,
MRI

##  1 Introduction

Support vector machines (SVMs) and other kernel based methods ushered in a new
wave of supervised learning methods based on the fundamental insight that
challenging tasks in low dimensions may become easier when the data is lifted
to higher dimensional spaces (Boser et al., 1992; Cortes and Vapnik, 1995;
Hofmann et al., 2008), as illustrated in Figure 1. These methods, however,
become prohibitive when dealing with massive datasets in high dimensional
feature spaces as their space complexity grows at least quadratically with the
number of data points (Bordes et al., 2005; Novikov et al., 2016). Further,
SVMs are also known to be sensitive to the specific choice of the kernel
parameters and as a consequence the decision boundaries learnt are prone to
over-fitting (Burges, 1998; Bordes et al., 2005). These factors have
discouraged their successful adoption to tasks involving large datasets
comprising high resolution images where deep learning based methods have shown
to fare better (Litjens et al., 2017; Liu et al., 2017).

Tensor networks, also known as tensor trains, offer a different and more
efficient framework to dealing with such high dimensional spaces.
Fundamentally, tensor networks are factorisations of higher order 222Order of
tensors are also referred to as ranks. In this work we adhere to using order.
tensors into networks of lower order tensors (Fannes et al., 1992; Oseledets,
2011; Bridgeman and Chubb, 2017). The number of parameters needed to specify
an order-$N$ tensor using tensor networks can be drastically reduced, from
exponential to polynomial dependence on $N$ (Perez-Garcia et al., 2006). This
massive reduction in number of parameters using tensor networks has been
predominantly applied to better understand quantum wave functions (Shi et al.,
2006). They have also seen applications in data compression (Cichocki et al.,
2016), and recently to better understand the expressive power of deep learning
models (Cohen et al., 2016; Glasser et al., 2019).

Figure 1: Data that is only non-linearly separable in lower dimensions can
become linearly separable in higher dimensions, illustrated here for a simple
case of 2-dimensional data that becomes linearly separable when lifted into
3-dimensions. This is the underlying principle in kernel based methods and
tensor networks which learn linear decision boundaries in high dimensional
spaces.

Recently, there has been an increasing interest in using tensor networks for
supervised learning, specifically focused on image classification tasks
(Stoudenmire and Schwab, 2016; Klus and Gelß, 2019; Efthymiou et al., 2019;
Sun et al., 2020). These methods rely on transforming two dimensional input
images (order-2 tensors) into one dimensional vectors (order-1 tensors) to
obtain linear decision boundaries in high dimensional spaces. Due to the
constraint of flattening to obtain vector inputs these methods are constrained
to work with images of small spatial resolution ($12$x$12$ px to $28$x$28$
px). Several improved flattening strategies have been attempted to maximize
the retained spatial correlation between pixels as the correlation declines
exponentially causing loss of information in larger images (Stoudenmire and
Schwab, 2016; Efthymiou et al., 2019; Cheng et al., 2019; Sagan, 2012). For
small enough images (like in MNIST 333http://yann.lecun.com/exdb/mnist/ or
Fashion MNIST 444https://github.com/zalandoresearch/fashion-mnist datasets)
there is some residual correlation in the flattened images which can be
exploited by lifting the vector input to higher dimensions using tensor
networks. This might not be effective in the case of images with higher
spatial resolution. Further, the information lost by flattening of images in
medical imaging tasks can be critical as the predicted decisions could be
dependent on the global structure of the image.

In this work, we present a tensor network based model adapted for classifying
two- and three-dimensional medical images which can be of higher spatial
resolution. The method presented here relies on lifting small image (slice or
volume) patches to higher dimensions, performing tensor network operations on
them to obtain intermediate representations which are then hierarchically
aggregated to predict the final classification decisions. These small, local
regions can be treated as being locally orderless drawing parallels to the
classical theory of locally orderless images in Koenderink and Van Doorn
(1999). As with locally orderless image analysis, we propose to extract useful
representations of small image regions in higher dimensions and aggregate them
in a hierarchical manner resulting in our locally orderless tensor network
(LoTeNet) model. The proposed LoTeNet model is used to learn decision
functions in high dimensional spaces in a supervised learning set-up and is
optimized end-to-end by backpropagating the error signal through the tensor
network. The LoTeNet model builds on an earlier work that adapted tensor
networks for supervised machine learning in Stoudenmire and Schwab (2016), and
also has similarities to the tensor network model in Efthymiou et al. (2019).
Specifically, as with both these models, the LoTeNet model uses the matrix
product state (MPS) tensor networks (Perez-Garcia et al., 2006) to approximate
linear decisions in a supervised learning set-up. The proposed modifications –
comprising the use of the hierarchical aggregation of patch-level
representations of image data using tensor networks – yield a large
computational advantage to LoTeNet making it amenable to be used in medical
image analysis.

The performance of the LoTeNet model is demonstrated using experiments on
classifying three medical imaging datasets: 2D histopathology images in PCam
dataset (Veeling et al., 2018), 2D thoracic computed tomography (CT) slices
from LIDC-IDRI dataset (Armato III et al., 2004) and 3D magnetric resonance
imaging (MRI) from the OASIS dataset (Marcus et al., 2007). These experiments
show that the LoTeNet model fares comparably to relevant state-of-the-art deep
learning methods requiring the tuning of a single model hyperparameter while
utilising only a fraction of the graphics processing unit (GPU) memory when
compared to their convolutional neural network (CNN) counterparts. The work
presented here is an extension of an earlier version of the method which
operated on two-dimensional data and is published as a conference publication
(Selvan and Dam, 2020).
The key contributions in this work are:

1\. A novel tensor network based model for classifying 2D and 3D medical image
data

2\. Extending supervised learning with tensor networks to images of higher
resolution

3\. Validation of the method on three public datasets and comparison to
relevant methods

4\. Demonstrating competitive performance using little GPU resources

5\. Conceptual connections of tensor networks with deep neural networks
In the remainder of this manuscript we present the basics of tensor notation,
tensor network fundamentals and formulate the supervised image classification
task in Section 2. The components of the proposed LoTeNet model and the final
model itself are presented in Section 3. The proposed model is discussed in
relation to existing literature in Section 4. Experiments on the three
datasets and results are presented in Section 5, along with the discussions
and future research directions in Section 6, and present overall conclusions
from this work in Section 7.

##  2 Background and Problem Formulation

We introduce key concepts pertaining the use and optimisation of tensor
networks in this section which will be put together to describe the proposed
method in the following sections.

Figure 2: Left: Tensor notation depicting a scalar $S$, vector ${V}^{i}$,
matrix $M^{ij}$ and a general order-$3$ tensor $T^{ijk}$. Center: Tensor
notation for matrix multiplication or _tensor contraction_ , which are used
extensively in the matrix product state networks used in this work. We adhere
to the convention that the contracted indices are written as subscripts.
Right: Tensor notation for trace of product of two matrices.

###  2.1 Tensor Network Notations

Tensor networks and operations on them are described using an intuitive
graphical notation, introduced in Penrose (1971). Figure 2 (left) shows the
commonly used notations for a scalar $S$, vector $V^{i}$, matrix $M^{ij}$ and
a general order-$3$ tensor $T^{ijk}$. The number of dimensions of a tensor is
captured by the number of edges emanating from the nodes denoted by the edge
indices. For instance, the vector $V^{i}$ is an order-$1$ tensor indicated by
the edge with index $i$ and an order-$3$ tensor has three indices $(i,j,k)$
depicted by the three edges, and so on.

Tensor notations can also be used to capture operations on higher-order
tensors succinctly as shown in Figure 2 (center) where matrix product is
depicted, which is also known as _tensor contraction_. The edge between the
tensor nodes $X^{i}_{j}$ and $Y^{k}_{j}$ is the dimension subsumed in matrix
multiplication resulting in the matrix $Z^{ik}$. More thorough introduction to
tensor notations can be found in Bridgeman and Chubb (2017).

###  2.2 Linear Model in Exponentially High Dimensions

A linear model in a sufficiently high dimensional space can be very powerful
(Novikov et al., 2016). In SVMs, this is accomplished by the _implicit_
mapping of the input data into an infinite dimensional space using radial
basis function kernels (Hofmann et al., 2008). In this section, we describe
the procedure followed in recent tensor network based works, including in
ours, to map the input data into an exponentially high dimensional space.

Consider an input vector ${\mathbf{x}}\in[0,1]^{N}$, which can be obtained by
flattening a 2D or 3D image with $N$ pixels with intensity values that are
normalized in the interval $[0,1]$. A commonly used feature map for tensor
networks is obtained from the _tensor product_ of pixel-wise feature maps
(Stoudenmire and Schwab, 2016):

$\Phi^{i_{1},i_{2},\dots
i_{N}}({\mathbf{x}})=\phi^{i_{1}}(x_{1})\otimes\phi^{i_{2}}(x_{2})\otimes\cdots\phi^{i_{N}}(x_{N})$
(1)

where the local feature map acting on a pixel $x_{j}$ is indicated by
$\phi^{i_{j}}(\cdot)$. The indices $i_{j}$ run from $1$ to $d$ where $d$ is
the dimension of the local feature map. The feature maps are usually simple
non-linear functions restricted to have unit norm such that the joint feature
map in Eq. (1) also has unit norm. A widely used local feature map with $d=2$
inspired from quantum wave function analysis is (Stoudenmire and Schwab, 2016)
is shown in (2), and a simpler local feature map from Efthymiou et al. (2019)
in Eq. (3) with similar properties (but non-unit norm):

$\displaystyle\phi^{i_{j}}(x_{j})$
$\displaystyle=[\cos\left(\frac{\pi}{2}x_{j}\right),\sin\left(\frac{\pi}{2}x_{j}\right)]$
(2) $\displaystyle\phi^{i_{j}}(x_{j})$ $\displaystyle=[x_{j},1-x_{j}].$ (3)

The joint feature map $\Phi({\mathbf{x}})$ in Eq. (1) is an order-$N$ tensor
due to tensor product of the $N$ order-$1$ local feature maps of dimension $d$
in Eq. (2). The joint feature map $\Phi({\mathbf{x}})$ can be seen as mapping
each image to a vector in the $d^{N}$ dimensional feature space (Stoudenmire
and Schwab, 2016). For multi-channel inputs, such as RGB images or other
imaging modalities, with $C$ input channels the local feature map can be
applied to each channel separately such that the resulting space is of
dimension $(d\cdot C)^{N}$ (Efthymiou et al., 2019).

Given this high dimensional joint feature map $\Phi({\mathbf{x}})$ of Eq. (1)
for the input data ${\mathbf{x}}$, a linear decision rule for a multi-class
classification task can be formulated as:

$\displaystyle f({\mathbf{x}})$
$\displaystyle=\arg\max_{m}f^{m}({\mathbf{x}}),$ (4) $\displaystyle
f^{m}({\mathbf{x}})$ $\displaystyle=W^{m}\cdot\Phi({\mathbf{x}}).$ (5)

where $m=[0,1,\dots M-1]$ are the $M$ classes, and the weight tensor $W^{m}$
is of order-($N+1$) with the output tensor index $m$.

The tensor notation representation of the linear model of Eq. (5) is shown in
Figure 3 (Step 1) where the first column of gray nodes are the individual
pixel feature maps of order-$1$ with feature dimension $d$. The local feature
maps are connected to the order-($N+1$) weight tensor $W^{m}$ along $N$ edges
with one edge for output of dimension $M$ marked with index $m$.

The order-($N$+1) weight tensor $W^{m}$ results in a total of $M\cdot d^{N}$
weight elements. Even for a relatively small gray scale image, say of size
$100\times 100$, the total number of weight elements in $W^{m}$ can be
massive: $2\cdot 2^{10000}\approx 10^{3010}$ which is exponentially more than
the number of atoms in the
universe555https://en.wikipedia.org/wiki/Observable_universe! The $d^{N}$ lift
in quantum physics is sometimes seen as a convenient illusion (Poulin et al.,
2011; Orús, 2014) as the most interesting behavior of systems can be captured
with fewer degrees of freedom in this high dimensional space. This is
analogous to using fewer than input feature dimensions after dimensionality
reduction operations. In the next section we will see how tensor networks can
access useful sub-spaces represented by such high dimensional tensors, with
parameters that grow linearly with $N$ instead of growing exponentially with
$N$.

###  2.3 Matrix Product State (MPS)

Figure 3: (Step 1) Linear model of Eq. (5) in tensor notation. (Step 2) MPS
approximation of the linear model. (Step 3) Series of tensor contractions done
with MPS to compute $W^{m}\cdot\Phi({\mathbf{x}})$ in Eq. (5)

Consider two vectors, $T^{i}$ and $U^{j}$ with indices $i$ and $j$
respectively. The tensor product 666Note that the tensor product of two
order-1 tensors is the same as the vector outer product. of these two vectors
yields an order-$2$ tensor or a matrix $X^{ij}$. The matrix product state
(MPS) (Perez-Garcia et al., 2006) is a type of tensor network that expands on
this notion of tensor products allowing the factorization of an order-$N$
tensor (with $N$ edges) into a chain of lower-order tensors. Specifically, the
MPS tensor network can factorise an order-$N$ tensor with a chain of order-$3$
(with three edges) except on the borders where they are of order-$2$, as shown
in Figure 3 (Step 2). Consider a tensor of order-$N$ with indices
$i_{i},i_{2},\dots i_{N}$, using MPS it can be approximated using lower-order
tensors as

$W^{m,i_{1},i_{2},\dots
i_{N}}=\sum_{\alpha_{1},\alpha_{2},\dots\alpha_{N}}A^{i_{1}}_{\alpha_{1}}A^{i_{2}}_{\alpha_{1}\alpha_{2}}A^{i_{3}}_{\alpha_{2}\alpha_{3}}\dots
A^{m,i_{j}}_{\alpha_{j}\alpha_{j+1}}\dots A^{i_{N}}_{\alpha_{N}},$ (6)

where $A^{i_{j}}$ are the lower-order tensors. The subscript indices
$\alpha_{j}$ are virtual indices that are contracted and are of dimension
$\beta$ referred as the _bond dimension_. The components of these intermediate
lower-order tensors $A^{i_{j}}$ form the tunable parameters of the MPS
approximation. Note that any $N$ dimensional tensor can be represented exactly
using an MPS tensor network if $\beta=d^{N/2}$. In most applications, however,
$\beta$ is fixed to a small value or allowed to adapt dynamically when the MPS
tensor network is used to approximate higher-order tensors (Perez-Garcia et
al., 2006; Miller, 2019).

Figure 4:  (a): Illustration of MPS factorisation of an order-5 tensor
$W^{43253}$ into five tensors of lower order (up to order-3 ) based on
Equation (6). The bond dimension in this factorisation is $\beta=2$ seen as
the subscript indices which are contracted. The tensor $W^{43253}$ has
$4\times 3\times 2\times 5\times 3=360$ parameters whereas the MPS
approximation requires $8+12+8+20+6=54$ parameters. (b): MPS approximation of
$W^{43253}$ in tensor notation.

One of the drawbacks of using MPS tensor networks is that they operate along
one dimension (as a chain). This is the primary reason that two dimensional
image data has to be first flattened to a vector when working with tensor
networks such as the MPS. Tensor networks that can work on arbitrary graphs,
which might be more suitable for image data like the projected entangled pair
states (PEPS) (Verstraete and Cirac, 2004), are not as well understood and do
not yet have efficient algorithms like the MPS.
Dot product with MPS: The decision function in Eq. (5) comprising the dot
product of the order-($N$+1) weight tensor $W^{m}$ and the joint feature map
$\Phi({\mathbf{x}})$ in Eq. (1) can be efficiently computed using the MPS
approximation in Eq. (6), depicted in Figure 3 (Step 2). The order in which
tensor contractions are performed can yield a computationally efficient
algorithm. The original MPS algorithm (Perez-Garcia et al., 2006) starts from
one of the ends, contracts a pair of tensors to obtain a new tensor which is
then contracted with the next tensor and this process is repeated until the
output tensor is reached. The cost of this algorithm scales with
{$N\cdot\beta^{3}\cdot d$} when compared to the cost without the MPS
approximation which scales with $d^{N}$. In this work, we use the MPS
implementation in Miller (2019) which performs parallel contraction of the
horizontal edges and then proceeds to contract them vertically as depicted in
Figure 3 (Step 3). The MPS approximation also reduces the number of tunable
parameters by an exponential factor, from {$M\cdot d^{N}$} to {$M\cdot d\cdot
N\cdot\beta^{2}$}.

##  3 Methods

The primary contribution in this work is a tensor network based image
classification model that can handle high resolution images of both two and
three spatial dimensions. Several recent tensor network models for supervised
image classification tasks operate on small planar images and flatten them
into vectors with different raveling strategies at the expense of global image
structure (Stoudenmire and Schwab, 2016; Han et al., 2018; Efthymiou et al.,
2019). In contrast to these methods, we only flatten small regions of the
images which can be assumed to be locally orderless (Koenderink and Van Doorn,
1999) and obtain expressive representations of local image regions using
tensor contractions in high dimensional spaces. We process these locally
orderless image regions in a hierarchical fashion using layers of MPS blocks
in the final model, which we call the locally orderless tensor network
(LoTeNet), shown in Figure 6. We next describe the proposed LoTeNet model in
detail.

###  3.1 Squeeze operation

Figure 5: Squeeze operation with stride $k=3$ which reshapes a $6\times
6\times 1$ image patch into a vector of size $4$ with feature dimension d=$9$.

The proposed LoTeNet model operates on small image regions and aggregates them
with MPS operations at multiple levels. These small image patches at each
level are created using the squeeze operation in two steps, as illustrated in
Figure 5 for an order-$2$ tensor i.e an image with two spatial dimensions.

The input images are converted into vectors with inflated feature dimensions
by applying kernels of stride $k$ along each of the spatial dimensions.
Consider an input image $X^{hij}$ with $S=3$ spatial dimensions (height:H,
width:V, depth:D), $N$ pixels and $C=d$ channels as feature dimension:
$X^{hij}\in{\mathbb{R}}^{H\times V\times D\times C}$. In the first step, input
image is reshaped into smaller patches controlled by the stride $k$:
$X^{hij}\in{\mathbb{R}}^{(H/k)\times(V/k)\times(D/k)\times d}$ such that the
feature dimension increases to $d=C\cdot k^{S}$. The stride of the kernel $k$
decides the extent of reduction in spatial dimensions and the corresponding
increase in feature dimension of the squeezed image.

In the second step, the reshaped image patches with inflated feature
dimensions are flattened from order-$S$ tensors to order-$1$ tensors,
resulting in $X^{h}\equiv{\mathbf{x}}\in{\mathbb{R}}^{(N/k^{S})}$ with
$d=C\cdot k^{S}$. This flattening operation provides spatial information in
the feature dimension to LoTeNet, thus retaining additional image level
structure when compared to flattening the entire image into a single vector
${\mathbf{x}}\in{\mathbb{R}}^{N}$ with $d=C$ channels (Efthymiou et al.,
2019). Further, the increase in the feature dimension $d=C\cdot k^{S}$ makes
the tensor network more expressive as it increases the dimensionality of the
feature space (Stoudenmire and Schwab, 2016).

The transformations in dimensions due to the squeeze operation parameterised
by the stride $k$, $\psi(\cdot;k)$, are summarised below:

$\psi(\cdot;k):\\{X^{hij}\in{\mathbb{R}}^{H\times V\times D\times
d},d=C\\}\longrightarrow\\{{\mathbf{x}}\in{\mathbb{R}}^{(N/k^{S})\times
d},d=C\cdot k^{S}\\}.$ (7)

The squeeze operation can also be treated as a local feature map due to the
increase in feature dimensions, similar to the more formal pixel level feature
maps mentioned in Eq. (2) and Eq. (3), operating on image patches instead of
individual pixels. Finally note that the squeeze kernel stride can be
different at each level $l$ of the model denoted $k_{l}$.

In the current formulation, the images are assumed to be rectangular in each
plane. As the main purpose of the squeeze operation is to increase the feature
dimensions by flattening small local neighbourhoods, this can be achieved by
reshaping any small, non-square region of the image into a vector of fixed
feature dimension. The simplest strategy to work with non-square (circular
images for instance) would be to pad the regions around the borders to make
them rectangular. Another approach could be to partition non-square images
into cells using Voronoi tesselations (Du et al., 1999) and squeezing these
cells into vectors.

###  3.2 Locally orderless tensor network (LoTeNet)

Figure 6: The proposed locally orderless tensor network (LoTeNet) with $L$
layers. Each layer $l$ consists of several MPS blocks based on the squeeze
kernel stride $k_{l}$ at that layer. The squeeze operation is as described in
Figure 5. The final MPS block outputs the $M$ class predictions shown as the
edge with index $m$.

An overview of the proposed LoTeNet model is shown in Figure 6 for an image
with two spatial dimensions i.e with $S=2$. The number of squeezed vectors and
the corresponding MPS blocks at each layer of the proposed model are also
indicated in Figure 6.

The LoTeNet model comprises of layers of MPS blocks interleaved with squeeze
operations without any non-linear components. The input to the first layer in
Figure 6 is the full resolution input image with $N$ pixels, shown with grids
marking the $k\times k$ patches which are then squeezed to obtain
$N_{1}=N/k^{S}$ vectors with $d=C\cdot k^{S}$. Each of these squeezed patches
are input into MPS blocks which contract the $d=C\cdot k^{S}$ vectors to
output a vector with dimension $d=\nu$. MPS blocks operating on squeezed image
patches can be interpreted as summarising them with a vector of size $\nu$
using a linear model in a higher dimensional feature space. In LoTeNet, we set
$\nu$ to be the same as the MPS bond dimension $\beta$, so that there is only
a single hyperparameter to be tuned.

The output vectors from all $N_{l}$ MPS blocks at a given layer $l$ are
reshaped back into the $S$ dimensional image space. However, due to the MPS
contractions, these intermediate image space representations will be of lower
resolution as indicated by the smaller image with fewer patches in Figure 6.
This is analogous to obtaining an average pooled version of the intermediate
feature maps in traditional CNNs. The intermediate images of lower resolution
formed at layer $l$ are further squeezed and contracted in the subsequent
layers of the model. This process is continued for $L-1$ layers and the final
MPS block acting as layer $L$ contracts the output vector from layer $L-1$ to
predict the decisions in Eq. (5). The weights of all the MPS blocks in each of
the $L$ layers are the model parameters which are tuned in a supervised
setting as described next.

###  3.3 Model Optimization

We view the sequence of MPS contractions in successive layers of our model as
forward propagation and rely on automatic differentiation to compute the
backward computation graph (Efthymiou et al., 2019). Torch MPS package
(Miller, 2019) is used to implement MPS blocks and trained in PyTorch (Paszke
et al., 2019) to learn the model parameters from training data in an end-to-
end fashion. These parameters are analogous to the weights of neural network
layers and can be updated in a similar iterative manner by backpropagating a
relevant error signal computed between the model predictions and the training
labels.

We minimize the cross-entropy loss between the true label
$y_{i}\in[0,\dots,M-1]$ for each image ${\mathbf{x}}_{i}\in\mathcal{D}$ and
the predicted label $f^{(y_{i})}({\mathbf{x}}_{i})$ in the training set
$\mathcal{D}$:

$\mathcal{L}(f^{(y_{i})}({\mathbf{x}}_{i}))=-\sum_{({\mathbf{x}}_{i},y_{i})\in\mathcal{D}}\log\frac{\exp{f^{(y_{i})}({\mathbf{x}}_{i})}}{\sum_{m=0}^{M-1}\exp{f^{(m)}({\mathbf{x}}_{i})}}=-\sum_{({\mathbf{x}}_{i},y_{i})\in\mathcal{D}}\log\Big{(}{\sigma}(f^{(y_{i})}({\mathbf{x}}_{i}))\Big{)}$
(8)

where $\sigma(\cdot)$ is the softmax operation used to obtain normalized
scores that can be interpreted as the predicted class probabilities. For
binary classes, the loss in Eq. (8) reduces to the binary cross entropy and
output dimension $M=1$ can be used with the sigmoid non-linearity to obtain
probabilistic predictions in $[0,1]$.

##  4 Related work

Key ideas in the proposed locally orderless tensor network are related to
several existing methods in the literature. The closest of them is the tensor
network model in Efthymiou et al. (2019) which in turn was based on the work
in Stoudenmire and Schwab (2016). The primary difference of LoTeNet compared
to Efthymiou et al. (2019) is in its ability to handle high resolution 2D/3D
images without losing spatial correlation between pixels. This capability is
incorporated with strategies from image analysis such as the interpretation of
small regions to be orderless (Koenderink and Van Doorn, 1999). Further, the
multi-layered approach used in LoTeNet has similarities with the multi-scale
analysis of images (Lindeberg, 2007).

The squeeze operation described in Section 3.1 serves dual purposes: to move
spatial information into feature dimension and to increase the feature
dimension. This is based on similar operations performed in normalizing flow
literature to provide additional spatial context via feature dimensions such
as in Dinh et al. (2017). The operation of stacking pixels into feature
dimension is also similar to the im2col operations used to transform
convolution into multiplications for improving efficiency of deep learning
operations (Chetlur et al., 2014). Recently, in Blendowski and Heinrich
(2020), the im2col operation was also used as a preprocessing step with
differentiable decision trees such as random ferns (Ozuysal et al., 2009). The
squeeze operation can also been as the inverse of the pixel shuffle operation
introduced in Shi et al. (2016) where the feature maps are interleaved into
spatial dimension to perform sub-pixel convolution.

##  5 Data and Experiments

###  5.1 Data

Figure 7: Two sample images from each of the three datasets a) PCam (Veeling
et al., 2018), b) LIDC (Armato III et al., 2004) and c) OASIS (Marcus et al.,
2007).

We report experiments on three public datasets for the task of binary
classification from 2D histopathology slides, 2D thoracic computed tomography
(CT) and 3D T1-weighted MRI scans.
PCam Dataset: The PatchCamelyon (PCam) dataset is a binary histopathology
image classification dataset introduced in Veeling et al. (2018). Image
patches of size $96\times 96$ px are extracted from the Camelyon16 Challenge
dataset (Bejnordi et al., 2017) with positive label indicating the presence of
at least one pixel of tumour tissue in the central $32\times 32$ px region and
a negative label indicating absence of tumour, as shown in Figure 7(a). In
this work, a modified PCam dataset from the Kaggle
challenge777https://www.kaggle.com/c/histopathologic-cancer-detection is used
that removes duplicate image patches included in the original dataset. This
dataset consists of about $220,000$ patches for training and an independent
test set of about $57,500$ patches is provided for evaluating the models. Data
augmentation is performed during training with $0.5$ probability comprising of
rotation, horizontal and vertical flips.
LIDC Dataset: The LIDC-IDRI datatset comprises of 1018 thoracic CT images with
lesions annotated by four radiologists (Armato III et al., 2004). We use the
$128\times 128$ px 2D slices from Knegt (2018) (shown in Figure 7 (b))
yielding a total of $15,096$ patches. The segmentation masks from the four
raters are transformed to binary labels indicating the presence (if two or
more radiologists marked a tumour) or absence of tumours (if less than two
raters marked tumours). These binary labels capture a majority voting among
the radiologists. The dataset is split into $60:20:20$ splits for training,
validation and a hold-out test set.
OASIS Dataset: The OASIS dataset consists of T1-weighted MRI scans of 416
subjects aged between 18 to 96 years (Marcus et al., 2007). The scans are
graded with clinical dementia rating (CDR) into four classes: non-demented=0,
very mild dementia=0.5, mild dementia=1.0, moderate dementia=2.0. We create an
Alzheimer’s disease classification dataset according to Wen et al. (2020) with
scans of all subjects above 60 years resulting in a dataset with 155 subjects.
Binary labels are extracted from the CDR scores: cognitively normal (CN) if
CDR=0 and Alzheimer’s disease (AD) if CDR$>$0, yielding 82:73 split between CN
and AD cases. The MRI scans are preprocessed using the extensive pipeline
described in Wen et al. (2020) comprising bias field correction, non-linear
registration and skull stripping using the ClincaDL
package888https://clinicadl.readthedocs.io/en/latest/. After preprocessing we
obtain volumes of size 128x128x128 voxels (Figure 7-c).

###  5.2 Experiments and Results

####  5.2.1 Implementation

The model is implemented in PyTorch (Paszke et al., 2019) based on the Torch
MPS package (Miller, 2019). The architecture of LoTeNet model is kept fixed
across all the datasets. The proposed LoTeNet model is evaluated with $L=4$
layers, squeeze kernel of size $k_{l}=2$ except for the input layer where
$k_{1}=8$ in order to reduce the number of MPS blocks in the first layer. The
most critical hyperparameter of LoTeNet is its bond dimension $\beta$; it was
set to $\beta=5$ obtained from the range $[2,3\dots 20]$ based on the
performance on the PCam validation set. We set the virtual dimension $\nu$ to
be the same as the bond dimension. Note that increasing the depth in LoTeNet
does not translate to obtaining more complex decisions (like in neural
networks) as LoTeNet is a linear model. Number of layers are controlled by the
kernel stride and primarily result in reduced computation cost (Selvan et al.,
2020). Finally, to keep the number of tunable hyperparameters lower, we also
do not optimize $\beta$ and $L$ jointly which would yield a model similar to
an adaptive MPS (Stoudenmire and Schwab, 2016; Miller, 2019)

We use the Adam optimizer (Kingma and Ba, 2014) with a learning rate of
$5\times 10^{-4}$. We use a batch size of $512$ for experiments on the PCam
and LIDC datasets, and a batch size of $4$ for the OASIS dataset. We
incorporate batch normalisation (Ioffe and Szegedy, 2015) after each layer in
all the models and we observed this resulted in faster and more robust
convergence. All models were assumed to have converged if there was no
improvement in validation accuracy over $10$ consecutive epochs and the model
with the best validation performance was used to predict on the test set. All
experiments were run on a single Tesla K80 GPU with 12 GB memory. The
development and training of all models (including baselines) in this work is
estimated to produce 14.1 kg of CO2eq, equivalent to 120.1 km travelled by car
as measured by Carbontracker 999https://github.com/lfwa/carbontracker/
(Anthony et al., 2020).

####  5.2.2 Experiments on 2D data

Table 1: Performance comparison on PCam dataset (left) and LIDC dataset
(right). For all models, we specify the GPU memory utilisation in gigabytes.
Best AUC across all models is shown in boldface.

PCam Models | GPU(GB) | AUC
---|---|---
Rotation Eq-CNN | $11.0$ | ${\bf 0.963}$
Densenet | $10.5$ | $0.962$
LoTeNet (ours) | $0.8$ | $0.943$
Tensor Net-X ($\beta=10$) | $5.2$ | $0.908$

LIDC Models | GPU(GB) | AUC
---|---|---
LoTeNet (ours) | $0.7$ | ${\bf 0.874}$
Tensor Net-X ($\beta=10)$ | $4.5$ | $0.847$
Densenet | $10.5$ | $0.829$
Tensor Net-X ($\beta=5)$ | $1.5$ | $0.823$

PCam models: Performance of the LoTeNet model is compared to the rotation
equivariant CNNs in (Veeling et al., 2018) which is also state-of-the-art for
the PCam dataset. Additionally, we compare to the single layer MPS model with
local feature map of the form in Eq. (3) with $\beta=10$ from Efthymiou et al.
(2019) reported here as Tensor Net-X and to Densenet (Huang et al., 2017) with
$4$ layers and a growth rate of $12$. We compare the binary classification
performance of the models using the area under the ROC curve (AUC) as it is
not sensitive to arbitrary decision thresholds. The performance metrics are
shown in Table 1 (left) along with the maximum GPU memory utilisation for each
of the models. We observe the AUC on the test set attained by LoTeNet model
($0.943$) is comparable to the Rotation Eq. CNN ($0.963$) and Densenet models
($0.962$) on the PCam dataset with a drastic reduction in the maximum GPU
memory utilisation; a mere 0.8 GB when compared to upto 11 GB for the Rotation
Eq.CNN 101010In Veeling et al. (2018) the authors used 4x12GB GPUs which was
confirmed to R.Selvan by B.Veeling. and 10.5 GB for Densenet. We also notice a
considerable improvement when compared to Tensor Net-X in the attained AUC
(0.908).
LIDC models: Similar to the PCam dataset, performance of LoTeNet model is
compared using AUC and maximum GPU memory utilisation, and are reported in
Table 1 (right). We compare to Densenet and Tensor Net-X with two bond
dimensions ($\beta=5,10$) to highlight the influence of the bond dimension.
LoTeNet model fares better than the compared models with an AUC of $0.874$
whereas the Densenet model achieves $0.829$ and Tensor Net-X achieves $0.847$
($\beta=10$) and $0.823$ ($\beta=5$). The GPU memory utlisation follows a
similar trend as with the PCam models with LoTeNet requiring only $0.7$ GB.

####  5.2.3 Experiments on 3D data

LoTeNet used in the 3D experiments is identical to the one used in 2D
experiments except for the local feature map. For the 3D experiments, we do
not apply the local feature map in Eq. (2) but instead rely on the increase in
feature dimension due to the squeeze operation in Eq. (7). This reduces the
number of parameters of the LoTeNet model and the risk of overfitting to the
limited data in the OASIS dataset.

The experimental set-up closely follows from Wen et al. (2020) and we use
5-fold cross validation to report the balanced accuracy on each of the test
folds. Balanced accuracy (BA) is binary accuracy normalised by the class skew
and helps when the classes are imbalanced. We reimplemented the subject-level
3D CNN used in Wen et al. (2020) as one of the compared methods. The subject-
level CNN comprises 5 layers of 3D CNN with an initial feature map of $8$, max
pooling and rectified linear units (ReLU) with doubling of feature maps after
each layer, and an additional three linear layers with ReLUs to output the
final predictions. We also report two other baseline models based on 3D CNN
and a multi-layered perceptron (MLP) model. The CNN baseline model is similar
to the subject-level CNN but has an initial feature map size of $32$. The MLP
baseline has four layers to match the LoTeNet model, and ReLU non-linear
activation function.

Table 2 summarises the performance of LoTeNet along with the compared methods
for two types of inputs. As LoTeNet is agnostic to the dimensionality of the
input image due to the squeeze operations, we demonstrate the performance on
volume data (3D) and also on 19,840 slices extracted along the z-axis (2D)
from the volume data. For the 2D experiments we ensure there is no data
leakage by creating folds at subject level. The balanced accuracy across five
test folds are shown and we see a clear improvement ($0.71\pm 0.09$) for the
3D case compared to subject-level CNN model ($0.67\pm 0.08$). LoTeNet model
has a lower score in the 2D case, indicating the model is able to extract more
global information from the volume data than from 2D slices. Additionally, we
report the number of parameters, maximum GPU utilisation and time per training
epoch for all the methods. To show the influence of the number of parameters,
we compare to a large MLP baseline model with 78M parameters and observe this
increase does not necessarily improve the baseline model’s performance.

Table 2: Performance comparison on OASIS dataset with comparing methods using 3D and 2D inputs. The performance is reported as balanced accuracy (BA) averaged over 5-fold cross validation. Number of parameters, maximum GPU utilization (GPU) and computation time per training epoch (t) for all methods are also reported. OASIS Models | Input | # Param. | GPU (GB) | t (s) | Average BA
---|---|---|---|---|---
LoTeNet (ours) | 3D | 52M | $2.1$ | 45.3 | ${\bf 0.71\pm 0.09}$
Subject-level CNN | 3D | 1M | $8.8$ | 8.7 | $0.67\pm 0.08$
CNN Baseline | 3D | 6.4M | $11.5$ | 12.8 | $0.64\pm 0.05$
MLP Baseline | 3D | 78M | $4.5$ | 4.1 | $0.63\pm 0.03$
Densenet | 2D | 0.2M | $10.5$ | 80.1 | $0.67\pm 0.04$
LoTeNet (ours) | 2D | 0.4M | $0.7$ | 81.3 | $0.65\pm 0.03$

##  6 Discussion

In this section, we focus on discussing high level trends observed with
experiments reported in Section 5.2, conceptual connections to other methods,
and point out possible directions for future work.

###  6.1 Comparison with tensor network models

The proposed LoTeNet model is related to the tensor network models used for
supervised learning tasks, such as the ones in Stoudenmire and Schwab (2016);
Efthymiou et al. (2019). In Stoudenmire and Schwab (2016), the parameters of
the tensor network are learned from data using density matrix renormalization
group (DMRG) algorithm (McCulloch, 2007), whereas in LoTeNet we use automatic
differentiation based on the implementation in Miller (2019) similar to
Efthymiou et al. (2019) to optimise the parameters of the model in Eq. (6).

The proposed LoTeNet model also improves upon the computation cost of
approximating linear decisions compared to other tensor network models such as
in Efthymiou et al. (2019). Consider the computation cost of approximating the
linear decision in Eq. (5) for an input image with $N$ pixels of $d$ features
with a bond dimension $\beta$. For the tensor network in Efthymiou et al.
(2019) it is $\mathcal{O}\left(N\cdot d\cdot\beta^{2}\right)$. For LoTeNet,
the computation cost reduces exponentially with the number of layers used:
$\mathcal{O}\left(\frac{N}{k^{2\cdot L}}\cdot L\cdot k^{2}\cdot
d\cdot\beta^{2}\right)$. The spatial resolution of the input image is reduced
with successive layers in LoTeNet (Figure 6) captured as the $k^{2\cdot L}$
term in the denominator. The cost of performing MPS operations on patches of
size $k\times k$ in $L$ layers is the additional $L\cdot k^{2}$ term in the
numerator. However, for $L\geq 2$ and $k>1$, the exponential reduction in
computation cost with $L$ can translate into considerable computational
advantage when processing high resolution medical images. Further, the
computation cost of LoTeNet can be reduced without noticeable degradation in
performance by sharing MPS in each layer across all patches (Selvan et al.,
2020).

###  6.2 On Performance

Figure 8: Influence of varying the bond dimensions, $\beta$ on the LIDC
dataset. Note the y-axis is between [0.8-1.0] for better visualisation.

The performance of LoTeNet across the three datasets is on par with, or
superior, to the other methods as reported in Tables 1 and 2. This robustness
across tasks is remarkable because the LoTeNet architecture is fixed across
the tasks. This behaviour could be attributed to the single model
hyperparameter – bond dimension, $\beta$ – which is fixed ($\beta=5$) across
all experiments. Recollect that the bond dimension controls the quality of MPS
approximations of the corresponding operations in high dimensional spaces as
described in Section 2.3. Consistent with the reporting in earlier works
(Efthymiou et al., 2019), we also find that the value of $\beta$ after a
specific value does not affect the model’s performance. This is illustrated in
Figure 8, where we show the sensitivity of the validation set performance on
LIDC dataset for different bond dimensions. This is a desirable feature as the
models can be expected to yield comparable performance when applied to
different datasets without entailing dataset specific cost of elaborate
architecture search. Additionally, this robustness could also be attributed to
more efficient computation of massive number of parameters in LoTeNet
resulting in generalised decision rules, as reported in Table 2 where we
notice that the number of parameters in LoTeNet are higher than other methods
(except MLP baseline).

In Tables 1 and 2, we report the GPU memory requirement for each of the
models. LoTeNet requires only a fraction of the GPU memory utilised by the
corresponding baseline methods, including the Densenet (Huang et al., 2017) or
Rotation Eq-CNN models (Veeling et al., 2018). One reason for this drastic
reduction in GPU memory utilisation is because tensor networks do not maintain
massive intermediate feature maps. This behaviour is similar to feedforward
neural networks and unlike CNNs which use a large chunk of GPU memory mainly
to store intermediate feature maps (Rhu et al., 2016). Further, as LoTeNet is
based on contracting input data into smaller tensors in a hierarchical manner
the memory consumption with successive contracted layers does not increase.
This can be an important feature in medical imaging applications as larger
images and larger batch sizes can be processed without compromising the
quality of the learned decision rules.

###  6.3 Choice of local feature maps

The local feature maps in Equations (2) and (3) are simple transformations
which have been used in the tensor network literature (Stoudenmire and Schwab,
2016; Efthymiou et al., 2019). More recently, other intensity based local
feature maps such as wavelet transforms have been attempted for 1D signal
classification in Reyes and Stoudenmire (2020). While we use the local map in
Equation (2) for 2D data and squeeze operation based local map for 3D data,
experiments on other possible local feature maps were also performed. For
instance, we experimented with jet based feature maps (Larsen et al., 2012)
that compute spatial gradients at multiple scales, and also an MLP based local
feature map acting on individual pixels. Both these local feature maps did not
yield any considerable performance improvement and to adhere to existing
tensor network literature we used the simple sinusoidal feature maps which
also do not have additional hyperparameters.

###  6.4 Comparison with Neural Networks

Figure 9: Conceptual comparison of a feedforward neural network such as an MLP
with the proposed LoTeNet model. Notice the functional correspondence between
linear and MPS layers, and non-linear activation functions and the squeeze
operations between the MLP and LoTeNet models, respectively. The input image
is flattened to a vector to input to the MLP model, whereas for the LoTeNet
model we squeeze small image regions before flattening as described in Section
3.1.

Feed-forward neural networks: The class of tensor network models adapted for
machine learning tasks (Stoudenmire and Schwab, 2016; Efthymiou et al., 2019)
are conceptually most similar to feed-forward neural networks such as the MLP,
as both classes of models operate on vector inputs (see Figure 9). The
difference, however, arises in the type of decision boundaries optimised by
each class of model as illustrated in Figure 1. MLP based feed-forward
networks obtain non-linear decision boundaries using several layers of neurons
and non-linear activation functions (such as sigmoid, ReLU etc.). Tensor
network models, including the proposed model, optimise linear decision
boundaries in exponentially high dimensional spaces without using any non-
linear activation functions.
Convolutional neural networks: LoTeNet uses MPS blocks on image patches and
aggregates these tensor representations in a hierarchical manner to learn the
decision rule. While this could appear similar to the use of strided
convolution kernels in CNNs, it is indeed closer to feed-forward neural
networks. The use of an MPS block per patch perhaps can be interpreted as a
form of strided convolution but only in how the local operation is performed
on each patch. The primary difference is in the weight sharing; the weights of
MPS blocks are not shared across the image. Each $k\times k$ image region has
an MPS block acting on it. This is indicated in Figure 6, where we report the
number of MPS blocks used per layer.

###  6.5 Limitations and Future Work

The computation time reported in Table 2 shows it to be higher for LoTeNet
($45.3$s) compared to the CNN ($12.8$s) or MLP ($4.1$s) baselines as efficient
implementations of tensor contractions are not natively supported in
frameworks such as PyTorch. There are ongoing efforts to further improve
efficient implementations which can accelerate tensor network operations
(Fishman et al., 2020; Novikov et al., 2020). Considering that LoTeNet usually
converges between 10-20 epochs, the increase in computation time with the
current implementation is not substantial as the model can be trained even on
3D data easily under an hour.

One possible drawback of having a single hyperparameter $\beta$ controlling
the MPS approximations is the lack of granular control of the model capacity.
In the 2D OASIS experiments (Table 2), we further investigated the lower
average balanced accuracy score by trying out different $\beta$ values. This
experiment revealed that LoTeNet was either under-fitting ($\beta<4$) or over-
fitting ($4\leq\beta\leq 8$) to the training set. This could be attributed to
the fact that the number of parameters grow quadratically with $\beta$ and
these discrete jumps make it harder to obtain models of optimal capacity for
some tasks.

Another concern with tensor network based models is their tendency to be over
parameterised as the number of parameters scale as {$d\cdot N\cdot\beta^{2}$}.
The linear dependency on the number of pixels, $N$, for volumetric data
results in a model with massive number of parameters at the outset which is
further aggravated by the quadratic dependency on $\beta$. While in our
experiments, this has not had an adverse effect on the performance, we have
observed the model is capable of over-fitting on small training data quite
easily. This can be a concern when dealing with datasets with few samples.

One way of handling these problems would be to introduce weight-sharing
between MPS blocks at each layer so that the number of parameters can be
reduced and controlled better. A possible future work is to investigate the
effect of sharing MPS blocks per layer and study if equivariance to
translation can be induced in tensor network based models (Cohen and Welling,
2016).

Finally, the issue of extending LoTeNet type models beyond classification or
regression tasks to more diverse medical imaging tasks such as segmentation
and registration are open research questions waiting to be explored further.
Early work on using tensor networks for 2D segmentation tasks was recently
reported in Selvan et al. (2021).

##  7 Conclusion

We have presented a tensor network model for classification of 2D and 3D
medical images, of higher spatial resolutions. Using LoTeNet we have shown
that aggregating MPS approximations on small image regions in a hierarchical
manner retains global structure of the image data. With experiments on three
datasets we have shown such models can yield competitive performance on a
variety of classification tasks. We have measured the GPU memory requirement
of the proposed model and shown it to be a fraction of the memory utilisation
of state-of-the-art CNN based models when trained to attain comparable
performance. Finally, we have demonstrated that a single model architecture
with $\beta=5$ (tuned on one of the datasets) can perform well on other
datasets also, which can be an appealing feature as it reduces dataset
specific architecture search.

Acknowledgments

Data used in this work were provided in part by OASIS: Cross-Sectional:
Principal Investigators: D. Marcus, R, Buckner, J, Csernansky J. Morris; P50
AG05681, P01 AG03991, P01 AG026276, R01 AG021910, P20 MH071616, U24 RR021382.
We thank Eric Kim for permitting us to modify and use Figure 1. The authors
also thank all reviewers and editors for their insightful feedback.

Ethical Standards

The work follows appropriate ethical standards in conducting research and
writing the manuscript. All data used in this project are from open
repositories.

Conflicts of Interest

We declare we do not have conflicts of interest.

## References

  * Anthony et al. [2020] Lasse F. Wolff Anthony, Benjamin Kanding, and Raghavendra Selvan.  Carbontracker: Tracking and predicting the carbon footprint of training deep learning models.  ICML Workshop on Challenges in Deploying and monitoring Machine Learning Systems, July 2020.  arXiv:2007.03051.
  * Armato III et al. [2004] Samuel G Armato III, Geoffrey McLennan, Michael F McNitt-Gray, Charles R Meyer, David Yankelevitz, Denise R Aberle, Claudia I Henschke, Eric A Hoffman, Ella A Kazerooni, Heber MacMahon, et al.  Lung image database consortium: developing a resource for the medical imaging research community.  _Radiology_ , 232(3):739–748, 2004.
  * Bejnordi et al. [2017] Babak Ehteshami Bejnordi, Mitko Veta, Paul Johannes Van Diest, Bram Van Ginneken, Nico Karssemeijer, Geert Litjens, Jeroen AWM Van Der Laak, Meyke Hermsen, Quirine F Manson, Maschenka Balkenhol, et al.  Diagnostic assessment of deep learning algorithms for detection of lymph node metastases in women with breast cancer.  _Jama_ , 318(22):2199–2210, 2017.
  * Blendowski and Heinrich [2020] Max Blendowski and Mattias P. Heinrich.  Learning to map between ferns with differentiable binary embedding networks.  In _Medical Imaging with Deep Learning_ , 2020.  URL https://openreview.net/forum?id=sdSOFa3Y09.
  * Bordes et al. [2005] Antoine Bordes, Seyda Ertekin, Jason Weston, and Léon Bottou.  Fast kernel classifiers with online and active learning.  _Journal of Machine Learning Research_ , 6(Sep):1579–1619, 2005.
  * Boser et al. [1992] Bernhard E Boser, Isabelle M Guyon, and Vladimir N Vapnik.  A training algorithm for optimal margin classifiers.  In _Proceedings of the fifth annual workshop on Computational learning theory_ , pages 144–152, 1992.
  * Bridgeman and Chubb [2017] Jacob C Bridgeman and Christopher T Chubb.  Hand-waving and interpretive dance: an introductory course on tensor networks.  _Journal of Physics A: Mathematical and Theoretical_ , 50(22):223001, 2017.
  * Burges [1998] Christopher JC Burges.  A tutorial on support vector machines for pattern recognition.  _Data mining and knowledge discovery_ , 2(2):121–167, 1998.
  * Cheng et al. [2019] Song Cheng, Lei Wang, Tao Xiang, and Pan Zhang.  Tree tensor networks for generative modeling.  _Physical Review B_ , 99(15):155131, 2019.
  * Chetlur et al. [2014] Sharan Chetlur, Cliff Woolley, Philippe Vandermersch, Jonathan Cohen, John Tran, Bryan Catanzaro, and Evan Shelhamer.  cudnn: Efficient primitives for deep learning.  _arXiv preprint arXiv:1410.0759_ , 2014.
  * Cichocki et al. [2016] Andrzej Cichocki, Namgil Lee, Ivan Oseledets, Anh-Huy Phan, Qibin Zhao, Danilo P Mandic, et al.  Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions.  _Foundations and Trends® in Machine Learning_ , 9(4-5):249–429, 2016.
  * Cohen et al. [2016] Nadav Cohen, Or Sharir, and Amnon Shashua.  On the expressive power of deep learning: A tensor analysis.  In _Conference on Learning Theory_ , pages 698–728, 2016.
  * Cohen and Welling [2016] Taco Cohen and Max Welling.  Group equivariant convolutional networks.  In _International conference on machine learning_ , pages 2990–2999, 2016.
  * Cortes and Vapnik [1995] Corinna Cortes and Vladimir Vapnik.  Support-vector networks.  _Machine learning_ , 20(3):273–297, 1995.
  * Dinh et al. [2017] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using Real NVP.  In _International Conference on Learning Representation_ , 2017.
  * Du et al. [1999] Qiang Du, Vance Faber, and Max Gunzburger.  Centroidal voronoi tessellations: Applications and algorithms.  _SIAM review_ , 41(4):637–676, 1999.
  * Efthymiou et al. [2019] Stavros Efthymiou, Jack Hidary, and Stefan Leichenauer.  Tensornetwork for Machine Learning.  _arXiv preprint arXiv:1906.06329_ , 2019.
  * Fannes et al. [1992] Mark Fannes, Bruno Nachtergaele, and Reinhard F Werner.  Finitely correlated states on quantum spin chains.  _Communications in mathematical physics_ , 144(3):443–490, 1992.
  * Fishman et al. [2020] Matthew Fishman, Steven R White, and E Miles Stoudenmire.  The itensor software library for tensor network calculations.  _arXiv preprint arXiv:2007.14822_ , 2020.
  * Glasser et al. [2019] Ivan Glasser, Ryan Sweke, Nicola Pancotti, Jens Eisert, and Ignacio Cirac.  Expressive power of tensor-network factorizations for probabilistic modeling.  In _Advances in Neural Information Processing Systems_ , pages 1496–1508, 2019.
  * Han et al. [2018] Zhao-Yu Han, Jun Wang, Heng Fan, Lei Wang, and Pan Zhang.  Unsupervised generative modeling using matrix product states.  _Physical Review X_ , 8(3):031012, 2018.
  * Hofmann et al. [2008] Thomas Hofmann, Bernhard Schölkopf, and Alexander J Smola.  Kernel methods in machine learning.  _The annals of statistics_ , pages 1171–1220, 2008.
  * Huang et al. [2017] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 4700–4708, 2017.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  _arXiv preprint arXiv:1502.03167_ , 2015.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Klus and Gelß [2019] Stefan Klus and Patrick Gelß.  Tensor-based algorithms for image classification.  _Algorithms_ , 12(11):240, 2019.
  * Knegt [2018] Stefan Knegt.  A Probabilistic U-Net for segmentation of ambiguous images implemented in PyTorch.  https://github.com/stefanknegt/Probabilistic-Unet-Pytorch, 2018\.
  * Koenderink and Van Doorn [1999] Jan J Koenderink and Andrea J Van Doorn.  The structure of locally orderless images.  _International Journal of Computer Vision_ , 31(2-3):159–168, 1999.
  * Larsen et al. [2012] Anders Boesen Lindbo Larsen, Sune Darkner, Anders Lindbjerg Dahl, and Kim Steenstrup Pedersen.  Jet-based local image descriptors.  In _European Conference on Computer Vision_ , pages 638–650. Springer, 2012.
  * Lindeberg [2007] Tony Lindeberg.  Scale-space.  _Wiley Encyclopedia of Computer Science and Engineering_ , pages 2495–2504, 2007.
  * Litjens et al. [2017] Geert Litjens, Thijs Kooi, Babak Ehteshami Bejnordi, Arnaud Arindra Adiyoso Setio, Francesco Ciompi, Mohsen Ghafoorian, Jeroen Awm Van Der Laak, Bram Van Ginneken, and Clara I Sánchez.  A survey on deep learning in medical image analysis.  _Medical image analysis_ , 42:60–88, 2017.
  * Liu et al. [2017] Peng Liu, Kim-Kwang Raymond Choo, Lizhe Wang, and Fang Huang.  Svm or deep learning? a comparative study on remote sensing image classification.  _Soft Computing_ , 21(23):7053–7065, 2017.
  * Marcus et al. [2007] Daniel S Marcus, Tracy H Wang, Jamie Parker, John G Csernansky, John C Morris, and Randy L Buckner.  Open access series of imaging studies (oasis): cross-sectional mri data in young, middle aged, nondemented, and demented older adults.  _Journal of cognitive neuroscience_ , 19(9):1498–1507, 2007.
  * McCulloch [2007] Ian P McCulloch.  From density-matrix renormalization group to matrix product states.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2007(10):P10014, 2007.
  * Miller [2019] Jacob Miller.  Torchmps.  https://github.com/jemisjoky/torchmps, 2019.
  * Novikov et al. [2016] Alexander Novikov, Mikhail Trofimov, and Ivan Oseledets.  Exponential machines.  _arXiv preprint arXiv:1605.03795_ , 2016.
  * Novikov et al. [2020] Alexander Novikov, Pavel Izmailov, Valentin Khrulkov, Michael Figurnov, and Ivan V Oseledets.  Tensor train decomposition on tensorflow (t3f).  _Journal of Machine Learning Research_ , 21(30):1–7, 2020.
  * Orús [2014] Román Orús.  A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  _Annals of Physics_ , 349:117–158, 2014.
  * Oseledets [2011] Ivan V Oseledets.  Tensor-train decomposition.  _SIAM Journal on Scientific Computing_ , 33(5):2295–2317, 2011.
  * Ozuysal et al. [2009] Mustafa Ozuysal, Michael Calonder, Vincent Lepetit, and Pascal Fua.  Fast keypoint recognition using random ferns.  _IEEE transactions on pattern analysis and machine intelligence_ , 32(3):448–461, 2009.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems_ , pages 8024–8035, 2019.
  * Penrose [1971] Roger Penrose.  Applications of negative dimensional tensors.  _Combinatorial mathematics and its applications_ , 1:221–244, 1971.
  * Perez-Garcia et al. [2006] David Perez-Garcia, Frank Verstraete, Michael M Wolf, and J Ignacio Cirac.  Matrix product state representations.  _arXiv preprint quant-ph/0608197_ , 2006.
  * Poulin et al. [2011] David Poulin, Angie Qarry, Rolando Somma, and Frank Verstraete.  Quantum simulation of time-dependent hamiltonians and the convenient illusion of hilbert space.  _Physical review letters_ , 106(17), 2011.
  * Reyes and Stoudenmire [2020] Justin Reyes and Miles Stoudenmire.  A multi-scale tensor network architecture for classification and regression.  _arXiv preprint arXiv:2001.08286_ , 2020.
  * Rhu et al. [2016] Minsoo Rhu, Natalia Gimelshein, Jason Clemons, Arslan Zulfiqar, and Stephen W Keckler.  vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design.  In _2016 49th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO)_ , pages 1–13. IEEE, 2016.
  * Sagan [2012] Hans Sagan.  _Space-filling curves_.  Springer Science & Business Media, 2012.
  * Selvan and Dam [2020] Raghavendra Selvan and Erik B Dam.  Tensor networks for medical image classification.  In _International Conference on Medical Imaging with Deep Learning – Full Paper Track_ , volume 121 of _Proceedings of Machine Learning Research_ , pages 721–732. PMLR, 06–08 Jul 2020.  URL http://proceedings.mlr.press/v121/selvan20a.html.
  * Selvan et al. [2020] Raghavendra Selvan, Silas Ørting, and Erik B Dam.  Multi-layered tensor networks for image classification.  In _First NeurIPS Workshop on Quantum Tensor Networks in Machine Learning_ , 2020.
  * Selvan et al. [2021] Raghavendra Selvan, Erik B Dam, and Jens Petersen.  Segmenting two-dimensional structures with strided tensor networks.  In _27th International Conference on Information Processing in Medical Imaging (IPMI-2021), Bornholm, Denmark._ , 2021.
  * Shi et al. [2016] Wenzhe Shi, Jose Caballero, Ferenc Huszár, Johannes Totz, Andrew P Aitken, Rob Bishop, Daniel Rueckert, and Zehan Wang.  Real-time single image and video super-resolution using an efficient sub-pixel convolutional neural network.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 1874–1883, 2016.
  * Shi et al. [2006] Y-Y Shi, L-M Duan, and Guifre Vidal.  Classical simulation of quantum many-body systems with a tree tensor network.  _Physical review a_ , 74(2):022320, 2006.
  * Stoudenmire and Schwab [2016] Edwin Stoudenmire and David J Schwab.  Supervised learning with tensor networks.  In _Advances in Neural Information Processing Systems_ , pages 4799–4807, 2016.
  * Sun et al. [2020] Zheng-Zhi Sun, Cheng Peng, Ding Liu, Shi-Ju Ran, and Gang Su.  Generative tensor network classification model for supervised machine learning.  _Physical Review B_ , 101(7):075135, 2020.
  * Veeling et al. [2018] Bastiaan S Veeling, Jasper Linmans, Jim Winkens, Taco Cohen, and Max Welling.  Rotation equivariant CNNs for digital pathology.  In _International Conference on Medical image computing and computer-assisted intervention_ , pages 210–218. Springer, 2018.
  * Verstraete and Cirac [2004] Frank Verstraete and J Ignacio Cirac.  Renormalization algorithms for quantum-many body systems in two and higher dimensions.  _arXiv preprint cond-mat/0407066_ , 2004.
  * Wen et al. [2020] Junhao Wen, Elina Thibeau-Sutre, Mauricio Diaz-Melo, Jorge Samper-González, Alexandre Routier, Simona Bottani, Didier Dormont, Stanley Durrleman, Ninon Burgos, Olivier Colliot, et al.  Convolutional Neural Networks for Classification of Alzheimer’s Disease: Overview and Reproducible Evaluation.  _Medical Image Analysis_ , page 101694, 2020.