# Finslerian Ressiner-Nordstrom spacetime

Xin Li 1<EMAIL_ADDRESS>1Department of Physics, Chongqing University,
Chongqing 401331, China

###### Abstract

We have obtained Finslerian Ressiner-Nordstrom solution where it is asymptotic
to a Finsler spacetime with constant flag curvature while
$r\rightarrow\infty$. The covariant derivative of modified Einstein tensor in
Finslerian gravitational field equation for this solution is conserved. The
symmetry of the special Finslerian Ressiner-Nordstrom spacetime, namely,
Finsler spacetime with constant flag curvature, has been investigated. It
admits four independent Killing vectors. The Finslerian Ressiner-Nordstrom
solution differs from Ressiner-Nordstrom metric only in two dimensional
subspace. And our solution requires that its two dimensional subspace have
constant flag curvature. We have obtained eigenfunction of Finslerian Laplace
operator of “Finslerian sphere”, namely, a special subspace with positive
constant flag curvature. The eigenfunction is of the form
$\bar{Y}_{l}^{m}=Y_{l}^{m}+\epsilon^{2}(C_{l+2}^{m}Y_{l+2}^{m}+C_{l-2}^{m}Y_{l-2}^{m})$
in powers of Finslerian parameter $\epsilon$, where $C_{l+2}^{m}$ and
$C_{l-2}^{m}$ are constant. However, the eigenvalue depends on both $l$ and
$m$. The eigenvalues corresponded to $Y_{1}^{0}$ remain the same with
Riemannian Laplace operator and the eigenvalues corresponded to $Y_{1}^{\pm
1}$ are different. This fact just reflect the symmetry of “Finslerian sphere”,
which admits a $z$-axis rotational symmetry and breaks other symmetry of
Riemannian sphere. The eigenfunction of Finslerian Laplace operator implies
that monopolar and dipolar terms of multipole expansion of gravitational
potential are unchanged and other multipole terms are changed.

##  I Introduction

A black hole is a specific region of spacetime that has strong enough gravity
such that even the light cannot escape from it. The black hole physics has
been discussed intensively by physicists. The researches on black hole physics
origins from the exact solution of Einstein’s gravitational field equation. In
four dimensional spacetime, Schwarzschild and Kerr solutions are exact
solutions of Einstein vacuum field equation, which correspond to the spherical
symmetry and axis symmetry respectively. And Reissner-Nordstrom spacetime is a
solution that corresponds to the gravitational field generated by a charged
and spherical symmetric gravitational source. General form of the solutions
are called Kerr-Newmann spacetime, which corresponds to the gravitational
field generated by a charged and axis symmetric gravitational source Kerr-
Newmann . Inspiring from Kerr-Newmann spacetime, the no hair theorem of the
black hole is found by physicists Gravitation . It states that the black hole
only has three properties, namely, mass, electric charge and spin. The above
solutions or spacetimes are asymptotically flat. Schwarzschild-de Sitter
spacetime is a solution with cosmological horizons, and it is asymptotic to de
Sitter spacetime Bousso . Based on the Kerr-Newmann spactime, Hawking and
Bekenstein et al. involved the concept of entropy for black hole Benkenstein
and constructed theory of black hole thermodynamics, then proposed four laws
of black hole thermodynamics Bardeen . Hawking radiation Hawking radiation ,
as a important complement of black hole thermodynamics, exhibits quantum
feature of black hole.

Three physical processes could lead to formations of black hole. One is the
gravitational collapse of a heavy starCelotti . The second is the
gravitational collapse of a primordial overdensity in the early universeCarr1
; Carr2 ; Carr3 ; Clesse ; Dolgov . The third is high-energy collisions
Giddings . The merge of black hole binary could generate gravitational wave,
such phenomena has been observed by the Advanced LIGO detectors, i.e.
GW150914, GW151226, GW170104 and GW170814 GW150914 ; GW151226 ; GW170104 ;
GW170814 .

It is interesting and important to search more exact solutions of
gravitational field equation in four dimensional spacetime. And the black
holes that corresponded to these solutions are expected to be tested by the
gravitational wave detectors in the near future. Finsler geometry Book by Bao
is a natural generation of Riemannian geometry. The basic feature of Finsler
geometry is that its length element does not have quadratic restriction.
Generally, Finslerian extension of a given Riemannian spacetime has lower
symmetry than the Riemannian spacetimeWang ; Finsler PF . A typical example of
Finsler spacetime, i.e., Randers spacetime Randers , breaks rotational
symmetry and induces parity violation. By this basic feature of Finsler
geometry, Finsler geometry is used to describe violation of Lorentz invariance
Girelli ; Gibbons ; Finsler LI ; Kostelecky  and study the anisotropy of our
universe Kouretsis ; Finsler alphae .

We have suggested an anisotropic inflation model in which the background space
is taken to be a Randers space. This anisotropic inflation model could account
for power asymmetry of the cosmic microwave background (CMB)Finsler inflation
. In Finsler geometry, there is no unique extension of Riemannian geometric
objects, such as the connections and curvatureBook by Bao . This feature also
involves in searching the gravitational field equation in Finsler spacetime Li
Berwald ; Miron ; Rutz ; Vacaru ; Pfeifer Gravity . We have proposed the
gravitational field equation in Finsler spacetime and found a non-Riemannian
exact solution Finsler BH . This solution or metric is none other than the
Schwarzschild metric except for the change from the Riemann sphere to
“Finslerian sphere”. And interior solution for the Finslerian Schwarzschild
metric also exists. We have proved that covariant derivative of the Finslerian
gravitational field equation for the metric is conserved. It is interesting to
search a general solution for our Finslerian gravitational field equation and
test if black hole corresponded to the solution possess the three properties,
namely, mass, electric charge and spin or not. The spherical harmonics are
eigenfunction of Laplace operator for the Riemann sphere, which play an
important role in modern physics. Since our Finslerian Schwarzschild solution
admits a “Finslerian sphere”, it is worth investigating eigenfunction of
Laplace operator for “Finslerian sphere”. It is expected that symmetry of
“Finslerian sphere” has direct influence on its eigenfunction.

This paper is organized as follows. In section II, we first give a brief
introduction to Ressiner-Nordstrom metric. Then, we present Finslerian
Ressiner-Nordstrom solution in Finsler spacetime. We discuss the symmetry of
Finslerian Ressiner-Nordstrom metric at the last of section II. In sec III, we
present Finslerian Laplace operator of “Finslerian sphere”, and give the
eigenfunction and corresponded eigenvalue of Finslerian Laplace operator.
Conclusions and remarks are given in section IV.

##  II Exact solution of gravitational field equation in Finsler spacetime

###  II.1 Brief introduction to Ressiner-Nordstrom metric

In general relativity, the Ressiner-Nordstrom spacetime is given as follows

$ds^{2}=-fdt^{2}+f^{-1}dr^{2}+r^{2}d\Omega^{2}_{k},$ (1)

where $f=k-\frac{2GM}{r}-br^{2}+\frac{4\pi GQ^{2}}{r^{2}}$, and the metric
$d\Omega^{2}_{k}$ denotes the 2 dimensional metric with constant sectional
curvature $k$. Usually, after a reparameterization, we can set $k$ to be
$1,0,-1$. The Ressiner-Nordstrom is a solution of Einstein field equation,
i.e.,

$Ric_{\mu\nu}-g_{\mu\nu}S/2=8\pi GT_{\mu\nu},$ (2)

where $Ric_{\mu\nu}$ is Ricci tensor and $S=g^{\mu\nu}Ric_{\mu\nu}$ is the
scalar curvature. The energy momentum tensor of the Ressiner-Nordstrom
spacetime is given as

$T_{\mu\nu}=T^{em}_{\mu\nu}+T^{c}_{\mu\nu}$ (3)

where
$T^{em}_{\mu\nu}=\frac{Q^{2}}{2r^{4}}diag\\{f,-f,r^{2}g^{\omega}_{ij}\\}$
($g^{\omega}_{ij}$ is the metric of $d\Omega^{2}_{k}$)denotes the energy
momentum tensor of electromagnetic field and
$T^{c}_{\mu\nu}=-3bg_{\mu\nu}/8\pi G$. The Ressiner-Nordstrom spacetime will
reduce to four dimensional spacetime with constant curvature, namely, the de
Sitter spacetime, if $Q=M=0$ and $k=1$.

###  II.2 Finslerian Ressiner-Nordstrom solution

Instead of defining an inner product structure over the tangent bundle in
Riemann geometry, Finsler geometry is based on the so called Finsler structure
$F$ with the property $F(x,\lambda y)=\lambda F(x,y)$ for all $\lambda>0$,
where $x\in M$ represents position and $y\equiv\frac{dx}{d\tau}$ represents
velocity. The Finslerian metric is given asBook by Bao

$g_{\mu\nu}\equiv\frac{\partial}{\partial y^{\mu}}\frac{\partial}{\partial
y^{\nu}}\left(\frac{1}{2}F^{2}\right).$ (4)

In our previous research Finsler BH , we have proposed that the vacuum
gravitational field equation is vanish of the Ricci scalar, and we have
obtained a solution of Finslerian vacuum field equation. It is given as

$F^{2}=-\left(1-\frac{2GM}{r}\right)y^{t}y^{t}+\left(1-\frac{2GM}{r}\right)^{-1}y^{r}y^{r}+r^{2}F^{2}_{FS},$
(5)

where $F_{FS}$ is a two dimensional Finsler spacetime with positive constant
flag curvature. The Ricci scalar is given as

$Ric\equiv R^{\mu}_{~{}\mu}=\frac{1}{F^{2}}\left(2\frac{\partial
G^{\mu}}{\partial x^{\mu}}-y^{\lambda}\frac{\partial^{2}G^{\mu}}{\partial
x^{\lambda}\partial y^{\mu}}+2G^{\lambda}\frac{\partial^{2}G^{\mu}}{\partial
y^{\lambda}\partial y^{\mu}}-\frac{\partial G^{\mu}}{\partial
y^{\lambda}}\frac{\partial G^{\lambda}}{\partial y^{\mu}}\right),$ (6)

where

$G^{\mu}=\frac{1}{4}g^{\mu\nu}\left(\frac{\partial^{2}F^{2}}{\partial
x^{\lambda}\partial y^{\nu}}y^{\lambda}-\frac{\partial F^{2}}{\partial
x^{\nu}}\right)$ (7)

is called geodesic spray coefficients.

Now, we propose an ansatz that the Finsler structure is of the form

$F^{2}=-f(r)y^{t}y^{t}+f(r)^{-1}y^{r}y^{r}+r^{2}\bar{F}^{2}(\theta,\varphi,y^{\theta},y^{\varphi}).$
(8)

Throughout the paper, the index labeled by Greek alphabet denote the index of
four dimensional spacetime $F$, and the index labeled by Latin alphabet denote
the index of two dimensional subspace $\bar{F}$. Plugging the Finsler
structure (8) into the formula (7), we obtain that

$\displaystyle G^{t}$ $\displaystyle=$
$\displaystyle\frac{f^{\prime}}{2f}y^{t}y^{r},$ (9) $\displaystyle G^{r}$
$\displaystyle=$
$\displaystyle-\frac{f^{\prime}}{4f}y^{r}y^{r}+\frac{ff^{\prime}}{4}y^{t}y^{t}-\frac{r}{2A}\bar{F}^{2},$
(10) $\displaystyle G^{\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{r}y^{\theta}y^{r}+\bar{G}^{\theta},$ (11)
$\displaystyle G^{\varphi}$ $\displaystyle=$
$\displaystyle\frac{1}{r}y^{\varphi}y^{r}+\bar{G}^{\varphi},$ (12)

where the prime denotes the derivative with respect to $r$, and the $\bar{G}$
is the geodesic spray coefficients derived by $\bar{F}$. Plugging the geodesic
coefficients (9,10,11,12) into the formula of Ricci scalar (6), we obtain that

$\displaystyle
F^{2}Ric=\left[\frac{ff^{\prime\prime}}{2}+\frac{ff^{\prime}}{r}\right]y^{t}y^{t}+\left[-\frac{f^{\prime\prime}}{2f}-\frac{f^{\prime}}{rf}\right]y^{r}y^{r}+\left[\bar{R}ic-
f-rf^{\prime}\right]\bar{F}^{2}~{}$ (13)

where $\bar{R}ic$ denotes the Ricci scalar of Finsler structure $\bar{F}$. We
have used the property of homogenous function $H(\lambda y)=\lambda^{n}H(y)$,
i.e. $y^{\mu}\frac{\partial H(y)}{\partial y^{\mu}}=nH(y)$, to derive the
geodesic spray coefficients $G^{\mu}$ and Ricci scalar. By equation (13), we
obtain the solution of $Ric=0$. It is of the form

$\displaystyle\bar{R}ic$ $\displaystyle=$ $\displaystyle k,$ (14)
$\displaystyle f_{S}$ $\displaystyle=$ $\displaystyle k-2GM/r,$ (15)

where $k=\pm 1,0$. And the solution of constant Ricci scalar ($Ric=3b$) is
given as

$\displaystyle\bar{R}ic$ $\displaystyle=$ $\displaystyle k,$ (16)
$\displaystyle f_{Sd}$ $\displaystyle=$ $\displaystyle k-2GM/r-br^{2}.$ (17)

A Finsler spacetime with constant flag curvature $K$ must have constant Ricci
scalar $(n-1)K$ ($n$ is the dimension of the spacetime). However, the reverse
statement is not true. Now, we test whether the solution (17) is corresponded
to the Finsler spacetime with constant flag curvature or not. A Finsler
spacetime with constant flag curvature $K$ is equivallent to its predecessor
of flag curvature possess the following form Book by Bao

$F^{2}R^{\mu}_{\nu}=K\left(F^{2}\delta^{\mu}_{\nu}-\frac{y^{\mu}}{2}\frac{\partial
F^{2}}{\partial y^{\nu}}\right),$ (18)

where $F^{2}R^{\mu}_{\nu}$ is defined as

$F^{2}R^{\mu}_{~{}\nu}=2\frac{\partial G^{\mu}}{\partial
x^{\nu}}-y^{\lambda}\frac{\partial^{2}G^{\mu}}{\partial x^{\lambda}\partial
y^{\nu}}+2G^{\lambda}\frac{\partial^{2}G^{\mu}}{\partial y^{\lambda}\partial
y^{\nu}}-\frac{\partial G^{\mu}}{\partial y^{\lambda}}\frac{\partial
G^{\lambda}}{\partial y^{\nu}}.$ (19)

Plugging the ansazt into the formula (19), after tedious calculation, we
obtain

$\displaystyle F^{2}R^{t}_{t}$ $\displaystyle=$
$\displaystyle-\frac{f^{\prime\prime}}{2f}y^{r}y^{r}-\frac{rf^{\prime}}{2}\bar{F}^{2},$
(20) $\displaystyle F^{2}R^{t}_{r}$ $\displaystyle=$
$\displaystyle\frac{f^{\prime\prime}}{2f}y^{t}y^{r},$ (21) $\displaystyle
F^{2}R^{t}_{i}$ $\displaystyle=$
$\displaystyle\frac{rf^{\prime}}{4}y^{t}\frac{\partial\bar{F}^{2}}{\partial
y^{i}},$ (22) $\displaystyle F^{2}R^{r}_{r}$ $\displaystyle=$
$\displaystyle\frac{ff^{\prime\prime}}{2}y^{t}y^{t}-\frac{rf^{\prime}}{2}\bar{F}^{2},$
(23) $\displaystyle F^{2}R^{r}_{i}$ $\displaystyle=$
$\displaystyle\frac{rf^{\prime}}{4}y^{r}\frac{\partial\bar{F}^{2}}{\partial
y^{i}},$ (24) $\displaystyle F^{2}R^{i}_{j}$ $\displaystyle=$
$\displaystyle\bar{F}^{2}\bar{R}^{i}_{j}+\left(\frac{ff^{\prime}}{2r}y^{t}y^{t}-\frac{f^{\prime}}{2rf}y^{r}y^{r}-f\bar{F}^{2}\right)\delta^{i}_{j}+\frac{f}{2}y^{i}\frac{\partial\bar{F}^{2}}{\partial
y^{j}}.$ (25)

Plugging the solution (17) into the above formula (20,21,22,23,24,25), it is
obvious from the formula (18) that Finsler metric (8) with the solution (17)
does not have constant flag curvature if the parameters $M\neq 0$ and $b\neq
0$. However, it is Einstein metric since its Ricci scalar is constant. By
making use of the property $R_{\mu\nu}=R_{\nu\mu}$ and noticing the $\bar{F}$
is a Finsler structure with constant Ricci scalar, we obtain that the special
case of the solution (17) with parameter $M=0$ and $b\neq 0$ is corresponded
to Finsler spacetime with constant flag curvature, i.e., it satisfies the
following relation

$F^{2}R^{\mu}_{\nu}=b\left(F^{2}\delta^{\mu}_{\nu}-\frac{y^{\mu}}{2}\frac{\partial
F^{2}}{\partial y^{\nu}}\right).$ (26)

It means that the flag curvature is constant and equals to $b$.

Pfeifer et al.have studied the electromagnetic field in Finsler spacetime
Pfeifer EM field . However, no specific solution of electromagnetic field
equation, such as static electric field, is discussed. Now, we study the
solution of Finslerian gravity with electric charge. Analogy to the Riemannian
Ressiner-Nordstrom metric, we suggest that the Finsler metric has the form of
the ansatz (8) with $f_{RN}=k-\frac{2GM}{r}+\frac{4\pi_{F}GQ^{2}}{r^{2}}$.
Then, by formula (13), its Ricci scalar is of the form

$\displaystyle
F^{2}Ric=\frac{4\pi_{F}GQ^{2}}{r^{4}}\left(fy^{t}y^{t}-f^{-1}y^{r}y^{r}+r^{2}\bar{F}^{2}\right).$
(27)

The gravitational field equation in the given Finsler spacetime (8) should be
of the form Finsler BH

$G^{\mu}_{\nu}=8\pi_{F}GT^{\mu}_{\nu},$ (28)

where the modified Einstein tensor $G_{\mu\nu}$ is defined as

$G_{\mu\nu}\equiv Ric_{\mu\nu}-\frac{1}{2}g_{\mu\nu}S,$ (29)

and $4\pi_{F}$ denotes the volume of $\bar{F}$. Here, the Ricci tensor we used
is first introduced by Akbar-ZadehAkbar

$Ric_{\mu\nu}=\frac{\partial^{2}\left(\frac{1}{2}F^{2}Ric\right)}{\partial
y^{\mu}\partial y^{\nu}}.$ (30)

And the scalar curvature in Finsler geometry is given as
$S=g^{\mu\nu}Ric_{\mu\nu}$. Then, by making use of the equations (27,29,30),
we obtain the non vanishing components of Einstein tensor

$Ric_{\mu\nu}=G_{\mu\nu}=\frac{4\pi_{F}GQ^{2}}{r^{4}}\rm
diag\\{f,-f^{-1},r^{2}\bar{g}_{ij}\\}.$ (31)

One should notice that the scalar curvature $S$ vanishes. It means that the
trace of energy momentum tensor vanishes. And such fact implies that the
particle is massless in physics. Plugging the result of Einstein tensor (31)
into the field equation (28), we obtain the energy momentum tensor

$T_{\mu\nu}=\frac{Q^{2}}{2r^{4}}\rm diag\\{f,-f^{-1},r^{2}\bar{g}_{ij}\\}.$
(32)

In ref. Finsler BH , we have proved that the covariant derivative of the
modified Einstein tensor is conserved in the Finsler spacetime (8), i.e.,
$G^{\mu}_{\nu|\mu}=0$, where “$|$” denotes the covariant derivative. Since the
Finslerian Schwarzschild-de Sitter solution $f=f_{Sd}$ and the Finslerian
Ressiner-Nordstrom solution $f=f_{RN}$ both possess the same form with
Finslerian Schwarzschild solution (8). Thus, following the same process of
ref. Finsler BH , one can find that the modified Einstein tensor and the
energy momentum tensor corresponded to the solution are also conserved. Energy
momentum tensor of general electromagnetic field in Finsler spacetime should
reduce to the energy momentum tensor given above (32) if the Finsler spacetime
reduce to Finslerian Ressiner-Nordstrom solution, i.e., the ansatz metric (8)
with $f=f_{RN}$.

###  II.3 Symmetry of Finslerian Ressiner-Nordstrom solution

In Riemannian geometry, spaces with constant sectional curvature are
equivalent. And all these spaces have $n(n+1)/2$ independent Killing vectors.
However, The number of independent Killing vectors of an $n$ dimensional non-
Riemannian Finsler spacetime should be no more than
$\frac{n(n-1)}{2}+1,~{}n\neq 4~{}\&~{}n\geq 3$ Wang . And a four dimensional
non-Riemannian Finsler spacetime has no more than 8 independent Killing
vectors Bogoslovsky . In general, Finslerian extension of a given Riemannian
spacetime has lower symmetry than the Riemannian spacetime. For example, the
Finslerian Schwarzschild spacetime $f=f_{S}$ only has two independent Killing
vectorsFinsler BH . The Kerr spacetime has two independent Killing vectors.
Thus, Finslerian extension of Kerr spacetime should only has one or zero
independent Killing vector. Furthermore, if we require the Finslerian
extensional spacetime is static, then, arbitrary static Finsler spacetime
should be Finslerian extension of Kerr spacetime. This fact implies that it is
hard to find Finslerian extension of Kerr spacetime. In Ref.Finsler BH , we
have shown that the Finsler spacetime with specific form (8) could form a
horizon at $f=0$. Therefore, at present, one can find from Finslerian
Ressiner-Nordstrom solution that Finslerian black hole has two properties,
namely, mass and electric charge.

The Killing equation $K_{V}(F)$ in Finsler spacetime is of the formFinsler PF

$K_{V}(F)\equiv\tilde{V}^{\mu}\frac{\partial F}{\partial
x^{\mu}}+y^{\nu}\frac{\partial\tilde{V}^{\mu}}{\partial x^{\nu}}\frac{\partial
F}{\partial y^{\mu}}=0.$ (33)

Plugging the formula (4) into the Killing equation (33), we obtain that

$V^{\mu}\frac{\partial g_{\alpha\beta}}{\partial
x^{\mu}}+g_{\alpha\lambda}\frac{\partial V^{\lambda}}{\partial
x^{\beta}}+g_{\lambda\beta}\frac{\partial V^{\lambda}}{\partial
x^{\alpha}}+y^{\nu}\frac{\partial V^{\mu}}{\partial x^{\nu}}\frac{\partial
g_{\alpha\beta}}{\partial y^{\mu}}=0.$ (34)

The left side of Killing equation (34) is just Lie derivative of Finsler
metric $g_{\alpha\beta}$ Yano book . From the equation (34), we investigate
the symmetry of the Finsler spacetime (8) with constant flag curvature, i.e.,
$f=f_{d}=k-br^{2}$. It should be noticed that the Killing equation (34)
differs from Riemannian Killing equation in $y^{\nu}\frac{\partial
V^{\mu}}{\partial x^{\nu}}\frac{\partial g_{\alpha\beta}}{\partial y^{\mu}}$.
And only component $g_{ij}$ of the Finsler metric of the Finsler spacetime (8)
has $y$ dependence. Also, Finsler spacetime (8) with $f=f_{d}$ reduces to
Riemannian spacetime with constant sectional curvature if $\bar{F}$ reduces to
Riemannian surface with constant sectional curvature. Therefore, there are
three independent Killing vectors which have index $t$ and $r$ only, and these
Killing vectors only depend on coordinate $t$ and $r$. And we have shown in
Ref.Finsler BH  that a “Finslerian sphere” with constant positive flag
curvature admits one Killing vector $V^{\varphi}=C^{\varphi}$, where
$C^{\varphi}$ is a constant. It is obvious that $V^{\varphi}=C^{\varphi}$ is
Killing vector of the Finsler spacetime (8) with constant flag curvature.
Finally, we conclude that the Finsler spacetime

$F^{2}_{d}=-(1-br^{2})y^{t}y^{t}+(1-br^{2})^{-1}y^{r}y^{r}+r^{2}F^{2}_{FS}$
(35)

admits four independent Killing vectors. The specific form of “Finslerian
sphere” will be given in next section. We have shown in Ref.Finsler PF  that a
four dimensional projectively flat Randers spacetime with constant flag
curvature admits six independent Killing vectors. Thus, Finsler spacetime (35)
and projectively flat Randers spacetime with constant flag curvature are not
equivalent, namely, after a coordinate transformation, one can change as
other. This fact is quite different with the one in Riemannian geometry. Since
each Riemannian spacetime with constant sectional curvature are equivalent.

##  III Finslerian Laplace operator on “Finslerian sphere”

Spherical harmonics is the eigenfunction of Laplace operator of Riemannian
2-sphere. Since our Finslerian Ressiner-Nordstrom solution in Finsler
spacetime differs from Ressiner-Nordstrom metric only in $\bar{F}$. It is
worth to investigate the Laplace operator of Finslerian surface with positive
constant flag curvature. In the discussion of above section, we have shown
that Finsler spaces with constant flag curvature may not equivalent to each
other. Thus, we adopt a specific form of Finsler surface with positive
constant flag curvature to study its Laplace operator, i.e., a two dimensional
Randers-Finsler space with constant positive flag curvature Bao

$F_{\rm
FS}=\frac{\sqrt{(1-\epsilon^{2}\sin^{2}\theta)y^{\theta}y^{\theta}+\sin^{2}\theta
y^{\varphi}y^{\varphi}}}{1-\epsilon^{2}\sin^{2}\theta}-\frac{\epsilon\sin^{2}\theta
y^{\varphi}}{1-\epsilon^{2}\sin^{2}\theta},$ (36)

where $0\leq\epsilon<1$. We call it “Finslerian sphere”.

In Riemannian geometry, the Laplace operator can be defined in several
different ways Gallot  and these Laplace operators are equivalent to each
other. However, Finslerian extension of these definitions will lead to
different Laplace operators. Various Finslerian Laplace operators are
respectively defined by Bao and LackeyBao and Lackey , ShenShen's book ,
Barthelme Barthelme . In this paper, we will adopt the Finslerian Laplace
operator defined by Barthelme. The Finslerian Laplace operator for “Finslerian
sphere” (36) is of the form Barthelme

$\displaystyle\Delta_{FS}$ $\displaystyle=$
$\displaystyle\frac{2(1-\epsilon^{2}\sin^{2}\theta)^{3/2}}{\sin^{2}\theta(1+\sqrt{1-\epsilon^{2}\sin^{2}\theta})}\frac{\partial^{2}}{\partial\varphi^{2}}+\frac{2(1-\epsilon^{2}\sin^{2}\theta)}{1+\sqrt{1-\epsilon^{2}\sin^{2}\theta}}\frac{\partial^{2}}{\partial\theta^{2}}$
(37)
$\displaystyle+\frac{2\cos\theta(\epsilon^{2}\sin^{2}\theta+\sqrt{1-\epsilon^{2}\sin^{2}\theta})}{\sin\theta(1+\sqrt{1-\epsilon^{2}\sin^{2}\theta})}\frac{\partial}{\partial\theta}.$

While $\epsilon=0$, the Finslerian Laplace operator (37) reduces to Riemanian
Laplace operator for 2-sphere. And its eigenfunction is just spherical
harmonics $Y_{lm}(\theta,\varphi)$ which satisfies

$\Delta_{FS}|_{\epsilon=0}Y_{l}^{m}=-l(l+1)Y_{l}^{m}.$ (38)

By using the equation (38), one can find that $Y_{0}^{0},Y_{1}^{0},Y_{1}^{\pm
1}$ are eigenfunctions of the Finslerian Laplace operator $\Delta_{FS}$, and
they satisfy the following equations

$\displaystyle\Delta_{FS}Y_{0}^{0}=0,~{}\Delta_{FS}Y_{1}^{0}=-2Y_{1}^{0},~{}\Delta_{FS}Y_{1}^{\pm
1}=(-2+2\epsilon^{2})Y_{1}^{\pm 1}~{}.$ (39)

It is obvious that the eigenvalues corresponded to $Y_{0}^{0},Y_{1}^{0}$
remain the same with Riemannian Laplace operator and the eigenvalues
corresponded to $Y_{1}^{\pm 1}$ are different. This fact just reflect the
symmetry of “Finslerian sphere”. As for the symmetry of “Finslerian sphere”,
We have showed that $z$-axis rotational symmetry, which corresponds to Killing
vector $V^{\varphi}=C^{\varphi}$, is preserved and other symmetry of Riemanian
2-sphere is broken Finsler BH .

The symmetry of “Finslerian sphere” may account for some specific physical
phenomena, such as the power asymmetry of CMB Power asymmetry . And such
phenomena can be treated as a perturbation of standard physical theories.
Thus, we expand Finslerian Laplace operator $\Delta_{FS}$ in powers of
$\epsilon$. To first order in $\epsilon^{2}$, the Finslerian Laplace operator
(37) is given as

$\displaystyle\Delta_{FS}$ $\displaystyle=$
$\displaystyle\frac{4-5\epsilon^{2}\sin^{2}\theta}{4\sin^{2}\theta}\frac{\partial^{2}}{\partial\varphi^{2}}+\left(1-\frac{3}{4}\epsilon^{2}\sin^{2}\theta\right)\frac{\partial^{2}}{\partial\theta^{2}}+\frac{\cos\theta}{\sin\theta}\left(1+\frac{3}{4}\epsilon^{2}\sin^{2}\theta\right)\frac{\partial}{\partial\theta}~{}.$
(40)

By making use of the recurrence formula of spherical harmonics, to first order
in $\epsilon^{2}$, we obtain the eigenfunction of the Finslerian Laplace
operator (40)

$\displaystyle\bar{Y}_{l}^{m}=Y_{l}^{m}+\epsilon^{2}(C_{l+2}^{m}Y_{l+2}^{m}+C_{l-2}^{m}Y_{l-2}^{m}),$
(41)

where

$\displaystyle C_{l+2}^{m}$ $\displaystyle=$
$\displaystyle-\frac{3l(l-1)}{8(2l+3)^{2}}\sqrt{\frac{(l+m+1)(l-m+1)(l+m+2)(l-m+2)}{(2l+1)(2l+5)}}~{},$
(42) $\displaystyle C_{l-2}^{m}$ $\displaystyle=$
$\displaystyle\frac{3(l+1)(l+2)}{8(2l-1)^{2}}\sqrt{\frac{(l+m)(l-m)(l+m-1)(l-m-1)}{(2l+1)(2l-3)}}~{}.$
(43)

And the corresponded eigenvalue of the Finslerian Laplace operator is given as

$\lambda=-l(l+1)+\epsilon^{2}\left(\frac{3(l-1)l(l+1)(l+2)}{2(2l-1)(2l+3)}+\frac{m^{2}(14l^{3}+21l^{2}+19l+6)}{2(2l+1)(2l-1)(2l+3)}\right)~{}.$
(44)

The spatial geometry of the Finslerian Schwarzschild metric is

$F^{2}_{3d}=y^{r}y^{r}+r^{2}\bar{F}^{2}_{FS},$ (45)

while $M=0$. Since $\bar{F}_{FS}$ does not depend on $r$, thus, following the
definition of Finslerian Laplace operator Barthelme , we find that Finslerian
Laplace equation for the three dimensional space (45) is of the form

$\Delta W=\frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial
W}{\partial r}\right)+\Delta_{FS}W=0~{}.$ (46)

Solution of the Finslerian Laplace equation (46) is of the form

$W=(Ar^{n_{1}}+Br^{n_{2}})\bar{Y}_{l}^{m},$ (47)

where $A$ and $B$ are constant that depend on boundary condition of the
Finslerian Laplace equation, and

$\displaystyle n_{1}$ $\displaystyle=$
$\displaystyle\frac{-1+\sqrt{1-4\lambda}}{2},$ (48) $\displaystyle n_{2}$
$\displaystyle=$ $\displaystyle\frac{-1-\sqrt{1-4\lambda}}{2}.$ (49)

Two facts for the solution (47) should be noticed. One is the index $n_{1}$
and $n_{2}$ depend not only on $l$ but also on $m$. Another fact is monopolar
and dipolar term of multipole expansion of gravitational potential is
unchanged in a Finsler spacetime with “Finslerian sphere”. It is consistent
with our previous research given in ref.Finsler BH  where we have shown that
the gravitational potential for the Schwarzschild-like spacetime is the same
with Newtonian gravity in weak field approximation. However, due to the
Finslerian modification of the eigenvalue $\lambda$ (44), other multipole
terms is changed.

##  IV Discussions and conclusions

Beside the definition of Ricci tensor introduced by Akbar-Zadeh, Shen
introduce another definition Shen Ric tensor

$Ric_{\mu\nu}=(R^{~{}\lambda}_{\mu~{}\lambda\nu}+R^{~{}\lambda}_{\nu~{}\lambda\mu})/2,$
(50)

where $R^{~{}\lambda}_{\mu~{}\lambda\nu}$ denotes the Riemann curvature tensor
of the Berwald connection. The two definitions are equivalent if the Finsler
spacetime has constant flag curvature. However, one can check that the two
definitions are not equivalent for the ansazt metric with
$f_{S}=k-\frac{2GM}{r}$, i.e., the Ricci flat case. Since the Ricci tensor
introduced by Akbar-Zadeh corresponded to Finslerian gravitational field
equation (28) with exact solutions, such as $f=f_{S}$, $f=f_{Sd}$ and
$f=f_{RN}$, therefore, the Ricci tensor introduced by Akbar-Zadeh is preferred
in physics.

In Finsler geometry, there are two types of volume form, namely, the Busemann-
Hausdorff volume form and Holmes-Thompson volume form Shen's book . We have
shown in Ref.Finsler BH  that the volume of “Finslerian sphere” in terms of
the Busemann-Hausdorff volume form is $4\pi$. The Finslerian Laplace operator
for “Finslerian sphere” (37) is defined on fiber bundle and its definition is
related to Holmes-Thompson volume form Barthelme . The volume of “Finslerian
sphere” in terms of Holmes-Thompson volume form is given as

${\rm
Vol}_{FS}=\int\frac{\sin\theta}{(1-\epsilon^{2}\sin^{2}\theta)^{3/2}}d\theta\wedge
d\varphi=\frac{4\pi}{1-\epsilon^{2}}.$ (51)

The two definitions of volume form will slightly alter the Finslerian
gravitation field equation (28), for the term $8\pi_{F}$ is double the surface
volume of Finsler space $\bar{F}$. In general relativity, black hole entropy
depends on the surface volume of the black hole. Thus, studying the black hole
thermodynamics in Finsler spacetime will help us to find which volume form is
preferred in physics. It will be discussed in our future work.

In this paper, we have obtained Finslerian Ressiner-Nordstrom solution where
it is asymptotic to a Finsler spacetime with constant flag curvature while
$r\rightarrow\infty$ (35). The covariant derivative of modified Einstein
tensor in Finslerian gravitational field equation for this solution is
conserved. The symmetry of the special Finslerian Ressiner-Nordstrom spacetime
(35) has been investigated. It admits four independent Killing vectors. The
Finslerian Ressiner-Nordstrom solution differs from Ressiner-Nordstrom metric
only in two dimensional subspace $\bar{F}$. The Finslerian Ressiner-Nordstrom
could form horizons at $f_{RN}=0$. At present, we can conclude that Finslerian
black hole have at least two properties, namely, mass and electric charge. Our
solutions show that two dimensional subspace $\bar{F}$ has constant flag
curvature. Spherical harmonics is the eigenfunction of Laplace operator of
Riemannian 2-sphere. We have obtained the eigenfunction of Finslerian Laplace
operator introduced by BarthelmeBarthelme  of “Finslerian sphere” (36). The
eigenfunction (41) is just a combination of spherical harmonics in powers of
Finslerian parameter $\epsilon$. However, the eigenvalue depends on both $l$
and $m$. The eigenvalues corresponded to $Y_{1}^{0}$ remain the same with
Riemannian Laplace operator and the eigenvalues corresponded to $Y_{1}^{\pm
1}$ are different. This fact just reflect the symmetry of “Finslerian sphere”,
which admits a Killing vector $V^{\varphi}=C^{\varphi}$ and breaks other
symmetry of Riemannian sphere. The eigenfunction (41) of Finslerian Laplace
operator implies that monopolar and dipolar terms of multipole expansion of
gravitational potential are unchanged and other multipole terms are changed.

###### Acknowledgements.

This work has been supported by the National Natural Science Fund of China
under grant Nos. 11775038.

## References

  * (1) E.T.Newman, A.I.Janis, J. Math. Phys. 6, 915 (1965).
  * (2) C.W. Misner, K.S. Thorne, and J.A. Wheeler, Gravitation, W.H. Freeman & Co., San Francisco, 1973.
  * (3) R. Bousso, arXiv:hep-th/0205177.
  * (4) J.D.Bekenstein, Phys.Rev.D 7, 2333 (1973).
  * (5) J. M. Bardeen, B. Carter, S.W. Hawking, Commun.Math.Phys. 31, 161 (1973).
  * (6) S.Hawking, Nature 248, 30 (1974); Commun.Math.Phys. 43, 199 (1975).
  * (7) A. Celotti, J.C. Miller, D.W. Sciama, Class. Quantum Gravity 16, A3 (1999).
  * (8) B. J. Carr and S. Hawking, Mon. Not. R. Astron. Soc. 168, 399 (1974).
  * (9) B. J. Carr, Astrophys. J. 201, 1 (1975).
  * (10) B. J. Carr, Astrophys. J. 206, 8 (1976).
  * (11) S. Clesse and J. Garca-Bellido, Phys. Rev. D 92, 023524 (2015).
  * (12) A. D. Dolgov and S. I. Blinnikov, Phys. Rev. D 89, 021301 (2014).
  * (13) Giddings, S. B.; Thomas, S. (2002),Physical Review D. 65 (5): 056010.
  * (14) B. P. Abbott et al (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016).
  * (15) B. P. Abbott et al (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 241103 (2016).
  * (16) B. P. Abbott et al (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 118, 221101 (2017).
  * (17) B. P. Abbott et al (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 141101 (2017).
  * (18) D. Bao, S. S. Chern, and Z. Shen, An Introduction to Riemann–Finsler Geometry, Graduate Texts in Mathematics 200, Springer, New York, 2000.
  * (19) H. C. Wang, J. Lond. Math. Soc. s1-22(1), 5 (1947).
  * (20) X. Li and Z. Chang, Differ. Geom. Appl. 30, 737 (2012).
  * (21) G. Randers, Phys. Rev. 59, 195 (1941).
  * (22) F. Girelli, S. Liberati, L. Sindoni, Phys. Rev. D 75, 064015 (2007).
  * (23) G.W. Gibbons, J. Gomis, and C. N. Pope, Phys. Rev. D 76, 081701 (2007).
  * (24) Z.Chang and X.Li, Phys. Lett. B 663, 103 (2008).
  * (25) V. A. Kostelecky, Phys. Lett. B 701, 137 (2011); V. A. Kostelecky, N. Russell, and R. Tsoc, Phys. Lett. B 716, 470 (2012).
  * (26) A. P. Kouretsis, M. Stathakopoulos, and P. C. Stavrinos, Phys. Rev. D 79, 104011 (2009); 82, 064035 (2010).
  * (27) Xin Li, Hai-Nan Lin, Sai Wang and Zhe Chang, Eur. Phys. J. C 75, 181 (2015).
  * (28) X. Li, S. Wang, Z. Chang, Eur. Phys. J. C 75, 260 (2015); X. Li, S. Wang, Eur. Phys. J. C 76, 51 (2016); X. Li, Hai-Nan Lin, Eur. Phys. J. C 77, 316 (2017).
  * (29) X. Li and Z. Chang, Chinese Physics C 34, 28 (2010).
  * (30) R. Miron and M. Anastasiei, The Geometry of Lagrange Spaces: Theory and Applications, Kluwer Acad. Publ. FTPH no. 59, (1994).
  * (31) S. F. Rutz, Gen. Relativ. Gravit. 25, 1139 (1993).
  * (32) S. Vacaru, P. Stavrinos, E. Gaburov, and D. Gonta, arXiv:gr-qc/0508023.
  * (33) C. Pfeifer and M. N. R. Wohlfarth, Phys. Rev. D 85, 064009 (2012).
  * (34) X. Li, Z. Chang, Phys. Rev. D 90, 064049 (2014).
  * (35) C. Pfeifer, M.N.R. Wohlfarth, Phys.Rev.D 84, 044039 (2011).
  * (36) H. Akbar-Zadeh, Acad. Roy. Belg. Bull. Cl. Sci. (5) 74, 281 (1988).
  * (37) G. Bogoslovsky, arXiv:0706.2621; G.W. Gibbons, J. Gomis, C.N. Pope, Phys. Rev. D 76, 081701 (2007).
  * (38) K. Yano, The theory of Lie derivative and its applications (North Holiland Amsterdam, 1955).
  * (39) D. Bao, C. Robles, and Z. Shen, J. Diff. Geom. 66, 377 (2004).
  * (40) S. Gallot, D. Hulin and J. Lafontaine, Riemannian Geometry, third edn., Universitext, Springer-Verlag, Berlin, 2004
  * (41) D.Bao, B. Lackey, C. R. Acad. Sc. Paris, 1996, 223(1): 51-56.
  * (42) Z. Shen, Lectures on Finsler Geometry, World Scientific, 2001.
  * (43) T. Barthelme, Israel J. Math., 2013, 196(1): 375-412
  * (44) Planck Collaboration, Astron. Astrophys. 571, A16 (2014).
  * (45) B. Li and Z. Shen, Chin.Phys.Lett. 20, 350 (2015).


© 2021 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

DOI:10.1109/ICIEA52957.2021.9436749

# Larger Sparse Quadratic Assignment Problem Optimization
Using Quantum Annealing and a Bit-Flip Heuristic Algorithm

Michiya Kuramata  Department of Industrial
Engineering and Economics
Tokyo Institute of Technology
Tokyo, Japan
Email<EMAIL_ADDRESS>Ryota Katsuki  NTT DATA Corporation
Tokyo, Japan
Email<EMAIL_ADDRESS>Kazuhide Nakata  Department of Industrial
Engineering and Economics
Tokyo Institute of Technology
Tokyo, Japan
Email<EMAIL_ADDRESS>
###### Abstract

Quantum annealing and D-Wave quantum annealer attracted considerable attention
for their ability to solve combinatorial optimization problems. In order to
solve other type of optimization problems, it is necessary to apply certain
kinds of mathematical transformations. However, linear constraints reduce the
size of problems that can be represented in quantum annealers, owing to the
sparseness of connections between qubits. For example, the quadratic
assignment problem (QAP) with linear equality constraints can be solved only
up to size 12 in the quantum annealer D-Wave Advantage, which has 5640 qubits.
To overcome this obstacle, Ohzeki developed a method for relaxing the linear
equality constraints and numerically verified the effectiveness of this method
with some target problems, but others remain unsolvable. In particular, it is
difficult to obtain feasible solutions to problems with hard constraints, such
as the QAP. We therefore propose a method for solving the QAP with quantum
annealing by applying a post-processing bit-flip heuristic algorithm to
solutions obtained by the Ohzeki method. In a numerical experiment, we solved
a sparse QAP by the proposed method. This sparse QAP has been used in areas
such as item listing on an E-commerce website. We successfully solved a QAP of
size 19 with high accuracy for the first time using D-Wave Advantage. We also
confirmed that the bit-flip heuristic algorithm moves infeasible solutions to
nearby feasible solutions in terms of Hamming distance with good computational
efficiency.

###### Index Terms:

optimization, quantum annealing, quadratic assignment problem

##  I Introduction

Quantum computers are expected to accelerate computation in fields such as
combinatorial optimization, simulation, machine learning, and so on. There are
two main types of quantum computers: gate-model quantum computers and quantum
annealers. In this study, we focus on quantum annealers that are hardware
implementations of quantum annealing [1]. Quantum annealers have started to
gain practical computational capabilities, particularly in combinatorial
optimization and sampling. D-Wave Systems has been developing quantum
annealers with steadily increasing numbers of qubits. In 2011, D-Wave Systems
released the first commercial quantum annealer, D-Wave One, which had 128
qubits. Since then, they released D-Wave Two with 512 qubits in 2013 and
D-Wave 2X with 1152 qubits in 2015, followed by the quantum annealer D-Wave
2000Q with 2048 qubits in 2017 and D-Wave Advantage [2] with 5640 qubits in
September 2020.

With this steady development of hardware, there has been an increasing
interest in finding use cases. Some practical applications to real-world
problems exist [3][4], particularly in the field of combinatorial
optimization. Neukart et al. [3] proposed an optimization problem for
assigning cars to roads to avoid congestion, which was solved using D-Wave 2X.
Nishimura et al. [4] used D-Wave 2000Q to optimize the order in which hotels
are listed on a Japanese hotel booking website (Jalan) as a quadratic
assignment problem (QAP) [5]. Both studies computed only small problems using
quantum annealers and larger problems using Qbsolv [6], an algorithm based on
tabu search.

The current D-Wave quantum annealer has many limitations. One is that there
are sparse connections between qubits, meaning there is no link between two
distant qubits, so multiple qubits may be required to represent a single
logical variable. Particularly in the case of a combinatorial optimization
problem such as QAP, which has a one-hot constraint such that only one of
multiple 0–1 binary variables takes a value of 1, a connection is required
between any pair of logical variables. The quantum annealer D-Wave Advantage
with 5640 qubits can thus only solve up to 180 fully connected logical
variables if missing qubits are not considered. Therefore, the maximum QAP
size that can be directly solved is 12 or less.

To overcome such obstacles, Ohzeki [7] proposed a method that relaxes linear
equality constraints such as one-hot constraints. Using this method, we can
represent larger problems on the quantum annealer hardware. Ohzeki numerically
verified the effectiveness of this method with some problems, but the target
problems cannot always be solved. In particular, it is difficult to obtain
feasible solutions of problems, such as QAP, as shown in Section IV.

##  II PREVIOUS STUDIES

###  II-A Quantum Annealing Using D-Wave Hardware

Like quantum annealing, simulated annealing (SA) [8] is a method for solving
combinatorial optimization problems. SA incorporates a pseudo-concept of heat,
which represents the ease of state change on a computer. Heat is initially
high to promote the active state transition. Then, the solution with the
minimum objective function value is found by gradually decreasing the heat. In
contrast, quantum annealing [1] obtains the minimum value of an objective
function by introducing and gradually weakening a transverse magnetic field
instead of SA heat. Quantum annealing is not performed pseudo-operatively on a
classical computer, but as a physical phenomenon in a quantum annealer, which
performs quantum annealing according to a physical model called the transverse
field Ising model, represented by

$\begin{split}\mathcal{H}(t)&=\left(1-s(t)\right)\left(-\sum_{i,j}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}-h\sum_{i}\sigma_{i}^{z}\right)\\\
&+s(t)\left(-\sum_{i}\sigma_{i}^{x}\right),\end{split}$ (1)

where $J_{ij}$ represents interaction of the $i$-th and $j$-th qubits,
$\sigma_{i}^{z}$ is the Pauli Z operator in the $i$-th qubit, $\sigma_{i}^{x}$
represents the Pauli X operator in the $i$-th qubit, and $s(t)$ is a function
of time $t$ that monotonically decreases from $1$ to $0$. In quantum
annealing, decreasing $s(t)$ corresponds to gradual weakening of the
transverse magnetic field, which finally gives the state $\\{s_{i}|i\in V\\}$
with a minimum value as an energy function of a physical model called the
Ising model, namely,

$E=\sum_{i,j}J_{ij}s_{i}s_{j}+\sum_{i}h_{i}s_{i},\quad s_{i}=\pm 1\;\;,\forall
i\in V,$ (2)

where $J_{ij}$ represents the coupling parameter of the $i$-th and $j$-th
qubits, $V$ is a set of qubits, and $h_{i}$ represents bias applied to the
$i$-th qubit. The combinatorial optimization problem can be solved using
physical hardware by setting $J_{ij}$ and $h_{i}$ so that the objective
function of the combinatorial optimization problem and the energy function of
the Ising model are equivalent. Handling combinatorial optimization problems
in the form of an Ising model is complicated because they are formulated in
terms of binary variables. Therefore, we represent the combinatorial
optimization problem as the quadratic unconstrained binary optimization (QUBO)
problem

$\begin{array}[]{ll}\text{minimize}&\mathbf{q}^{\top}Q\mathbf{q}\\\
\text{subject to}&\mathbf{q}\in\\{0,1\\}^{L},\end{array}$ (3)

where $q_{i}$ is a binary variable, $L$ is the number of $q_{i}$, and $Q$ is
an $L\times L$ matrix. QUBO is the problem of finding $\mathbf{q}$ that
minimizes $\mathbf{q}^{\top}Q\mathbf{q}$, which is equivalent to the Ising
model. In other words, a combinatorial optimization problem formulated as QUBO
can be optimized using quantum annealing. We can obtain QUBO, substituting
$s_{i}=2q_{i}-1$ into the energy function of the Ising model of (2). QUBO does
not allow for constraints. Therefore, to ensure that a solution satisfies its
constraints, QUBO needs a penalty term in the objective function. In this
paper, the combinatorial optimization problem is also represented as QUBO and
optimized by quantum annealing.

D-Wave Systems’ quantum processing unit (QPU) has a graph structure with
qubits as nodes and connections between qubits as edges. D-Wave Advantage uses
the Pegasus graph shown in Figure 1, where blue circles represent qubits and
lines represent connections between qubits. One qubit does not have a
connection with all other qubits, and connections between qubits are sparse.
Therefore, representing a single logical variable may require multiple qubits.
In other words, if we map multiple well-connected logical variables to the
Pegasus graph, we find that they consume more qubits to represent a single
logical variable. The number of connections between logical variables
corresponds to the number of nonzero elements of $Q$ in (3). Therefore, QUBO
with $Q$ with many nonzero elements is harder to represent in the Pegasus
graph, while QUBO with $Q$ with few nonzero elements is easier to represent.
Particularly in the case of a combinatorial optimization problem with one-hot
constraints, connections are needed between any pair of logical variables in
those constraints. Therefore, D-Wave Advantage cannot solve problems of size
13 or larger, such as the QAP or the traveling salesman problem (TSP).

Figure 1: Pegasus graph. Blue circles represent qubits, and lines connecting
qubits represent connections between them.

Note that in this paper, we define a one-hot constraint as one where exactly
one of multiple binary variables takes value 1. In a QAP, for example, the
constraint that a facility can be built at exactly one location is a one-hot
constraint. A feasible solution is one that satisfies all constraints, and an
infeasible solution is one that violates a constraint.

###  II-B The Quadratic Assignment Problem in E-commerce Websites

Nishimura et al. [4] proposed the following QAP for optimizing the order of
items on an E-commerce website, in consideration of their diversity:

$\begin{array}[]{ll}\text{ maximize
}&\sum_{i\in\text{Item}}\sum_{j\in\text{Site}}s_{ij}q_{ij}-w\\\
&\sum_{i\in\text{Item}}\sum_{i^{\prime}\in\text{Item}}\\\
&\sum_{j\in\text{Site}}\sum_{j^{\prime}\in\text{Site}}f_{ii^{\prime}}d_{jj^{\prime}}q_{ij}q_{i^{\prime}j^{\prime}}\\\
\text{ subject to }&\sum_{i\in\text{Item}}q_{ij}=1,\quad\forall
j\in\text{Site},\\\ &\sum_{j\in\text{Site}}q_{ij}=1,\quad\forall
i\in\text{Item},\\\ &q_{ij}\in\\{0,1\\},\quad\forall
i\in\text{Item},\quad\forall j\in\text{Site}.\end{array}$ (4)

Here the first term of the objective function represents estimated sales
without considering similarity, the second is a penalty term to avoid
sequences of highly similar items, “Item” is the set of items handled by the
website, and “Site” is a set of positions where items are placed on the
website. $q_{ij}$ is a binary variable that takes a value of 1 when item $i$
is placed in position $j$, $s_{ij}$ is estimated sales when item $i$ is in
position $j$, $f_{ii^{\prime}}$ is similarity between items $i$ and
$i^{\prime}$, $d_{jj^{\prime}}$ takes a value of 1 if positions $j$ and
$j^{\prime}$ are adjacent and 0 otherwise, and $w$ is a penalty parameter.

In other words, the first term of the objective function in (4) represents
estimated sales, and the second term penalizes orders of highly similar items.
The first constraint is that each item must be placed in one position, and the
second is that an item must occupy each position. By solving the optimization
problem in (4), it is possible to find an item arrangement that maximizes
sales while avoiding sequences of similar items. Formulating (4) as QUBO, we
obtain

$\begin{array}[]{ll}\text{minimize}&-\sum_{i\in\text{Item}}\sum_{j\in\text{Site}}s_{ij}q_{ij}\\\
&+w\sum_{i\in\text{Item}}\sum_{i^{\prime}\in\text{Item}}\\\
&\sum_{j\in\text{Site}}\sum_{j^{\prime}\in\text{Site}}f_{ii^{\prime}}d_{jj^{\prime}}q_{ij}q_{i^{\prime}j^{\prime}}\\\
&+\lambda\sum_{i\in\text{Item}}\left(\sum_{j\in\text{Site}}q_{ij}-1\right)^{2}\\\
&+\lambda\sum_{j\in\text{Site}}\left(\sum_{i\in\text{Item}}q_{ij}-1\right)^{2}\\\
\text{ subject to }&q_{ij}\in\\{0,1\\},\quad\forall
i\in\text{Item},\quad\forall j\in\text{Site},\end{array}$ (5)

where $\lambda$ is a penalty parameter with a large value that makes the
solution satisfy the constraints. We can optimize (5) by quantum annealing,
but as explained in Section II-A, the one-hot constraints make couplers
between logical variables fully connected, reducing the maximum size of
solvable problems.

###  II-C Breaking One-hot Constraints

In the QUBO formulation, we express constraints as a penalty term in the
objective function. Especially in problems with multiple one-hot constraints,
such as QAP and TSP, penalty terms create connections between arbitrary
logical variables. As described in Section II-A, the graph structure of a
quantum annealer has sparse connections between qubits, so it is necessary to
use multiple qubits to represent one logical variable. As a result, D-Wave
Advantage can only solve QAP up to size 12. This is a significant obstacle to
optimization using quantum annealing. Ohzeki [7] proposed a method, which we
call the Ohzeki method in this paper, for efficiently solving combinatorial
optimization problems with linear equation constraints such as one-hot
constraints using a quantum annealer. The Ohzeki method converts the quadratic
term generated by the penalty term to a linear term. As a result, the
connections between qubits become sparse, making it possible to solve larger
optimization problems with multiple one-hot constraints, such as QAP, with a
quantum annealer, but conversely making it difficult to obtain a feasible
solution. Briefly,

$\begin{array}[]{ll}\text{ minimize }&f_{0}(\mathbf{q})\\\ \text{ subject to
}&F_{k}(\mathbf{q})=C_{k},\quad\forall k\in\\{1,2,\ldots,M\\},\\\
&\mathbf{q}\in\\{0,1\\}^{L},\end{array}$ (6)

is a combinatorial optimization problem with equality constraints, where
$\mathbf{q}=(q_{1},q_{2},\ldots,q_{L})$ are binary variables, $L$ is the
number of binary variables, $M$ is the number of constraints,
$F_{k}(\mathbf{q})$ is a linear function of $\mathbf{q}$, and $f_{0}$ is the
objective function. The constraint $F_{k}(\mathbf{q})=C_{k}$ is a linear
equality constraint. We can regard this as one-hot constraint if we set
$F_{k}(\mathbf{q})=\sum_{k=1}^{L}q_{k}$, $C_{k}=1$. For optimization by a
quantum annealer, we transform (6) into the QUBO problem

$\begin{array}[]{ll}\text{minimize}&f(\mathbf{q})=f_{0}(\mathbf{q})+\frac{1}{2}\sum_{k=1}^{M}\lambda\left(F_{k}(\mathbf{q})-C_{k}\right)^{2}\\\
\text{ subject to }&\mathbf{q}\in\\{0,1\\}^{L},\end{array}$ (7)

where $\lambda$ is a penalty parameter, the second term represents a penalty
to satisfy the constraints, and $F_{k}(\mathbf{q})=C_{k}\;\;\forall k$.

The solution from a quantum annealer follows a Boltzmann distribution. We
applied the Hubbard–Stratonovich transformation [9][10] to the second term
$\frac{1}{2}\sum_{k=1}^{M}\lambda\left(F_{k}(\mathbf{q})-C_{k}\right)^{2}$ in
the objective function in (7), thus obtaining

$f_{\text{new}}(\mathbf{q};\mathbf{v})=f_{0}(\mathbf{q})-\sum_{k=1}^{M}v_{k}^{t}F_{k}(\mathbf{q}),$
(8)

where $v_{k}^{t}$ is a parameter and $t$ is an index representing the number
of iterations of the algorithm. This eliminates the dense connections between
logical variables caused by $\left(F_{k}(\mathbf{q})-C_{k}\right)^{2}$,
resolving the obstacle of reducing the problem size caused by the sparse QPU
graph structure (Figure 1).

The solution sampled by the quantum annealer for $f_{\text{new}}$ in (8)
follows the Boltzmann distribution

$Q(\mathbf{q})^{t}=\frac{1}{Z(\mathbf{v})}\exp{\left(-\beta
f_{\text{new}}\left(\mathbf{q};\mathbf{v}^{t}\right)\right)},$ (9)

with partition function

$Z(\mathbf{v})^{t}=\sum_{\mathbf{q}}\exp{\left(-\beta
f_{\text{new}}\left(\mathbf{q};\mathbf{v}^{t}\right)\right)},$ (10)

where $\beta$ is a reverse temperature parameter. Index $t$ in
$Q(\mathbf{q})^{t}$, $Z(\mathbf{v})^{t}$, and $\mathbf{v}^{t}$ represents
iteration $t$ of the algorithm. While this method eliminates the dense
connections between logical variables caused by the linear equality
constraints, doing so does not satisfy the linear equality constraint
$F_{k}(\mathbf{q})=C_{k}$ of the original combinatorial optimization problem
(6). To satisfy the constraints, we adjust $\mathbf{v}$ as

$v_{k}^{t+1}\leftarrow v_{k}^{t}+\eta\left(C_{k}-\left\langle
F_{k}(\mathbf{q})\right\rangle_{Q^{t}}\right),\quad\forall
k\in\\{1,2,\ldots,M\\},$ (11)

where $\left\langle F_{k}(\mathbf{q})\right\rangle_{Q^{t}}$ represents the
expected value of $F_{k}$ under distribution $Q^{t}$, and $\eta$ is a positive
parameter. Calculation of $\left\langle
F_{k}(\mathbf{q})\right\rangle_{Q^{t}}$ requires the expected value and is not
easily computed on a classical computer. However, we can use quantum annealers
to efficiently perform this computation. Since solutions
$(\mathbf{q}_{1},\mathbf{q}_{2},\ldots,\mathbf{q}_{N})$ are sampled from the
quantum annealer, $\left\langle F_{k}(\mathbf{q})\right\rangle_{Q^{t}}$ equals
the average $\frac{1}{N}\sum_{i=1}^{N}F_{k}(\mathbf{q}_{i})$, where $N$ is the
number of solutions sampled from the quantum annealer. We adjust $\mathbf{v}$
to produce feasible solutions according to (11), but as the experiments
described in Section IV show, all sampled solutions generated by the Ohzeki
method are infeasible despite satisfying some one-hot constraints. We
therefore propose a heuristic post-processing method in Section III-B.

##  III PROPOSED METHOD

###  III-A Applying the Ohzeki Method to the QAP on an E-commerce Website

We solve the optimization problem (4) on the E-commerce website described in
Section II-B with a size larger than that used in [4]. However, as mentioned
above, the optimization problem (4) is the QAP, and the one-hot constraint
causes the problem size to shrink. Therefore, in [4], problems with up to only
8 items can be solved by the quantum annealer D-Wave 2000Q. Even if D-Wave
Advantage is used in the same way, problems with 13 or more items cannot be
solved. To address this problem, we apply the Ohzeki method to (4). As a
result, while we did not obtain a feasible solution, we can embed a 20-item
problem (4) into D-Wave Advantage. Further, if using our proposed method as
described in Section III-B with the Ohzeki method, we can obtain approximate
solutions for problems with up to 19 items.

In this section, we describe how to apply the Ohzeki method to the
optimization problem (4) on an E-commerce website. First, we transform (5),
which is the QUBO representation of (4), into the notation of (7). We
respectively define $f_{0}$, $F_{i}^{\text{Item}}$, and $F_{j}^{\text{Site}}$
as (12), (13), and (13). $F_{i}^{\text{Item}}$ represents constraints on items
on the E-commerce website. Similarly, $F_{j}^{\text{Site}}$ represents
constraints on the position on the website. $F_{i}^{\text{Item}}$ and
$F_{j}^{\text{Site}}$ correspond to $F_{k}$ in (7), and $\mathbf{q}$ is an
array of binary variables. We define $n$ as the number of items and
$n=|\text{Item}|=|\text{Site}|$ holds. Using

$\begin{split}f_{0}(\mathbf{q})&=-\sum_{i\in\text{Item}}\sum_{j\in\text{Site}}s_{ij}q_{ij}\\\
&+w\sum_{i\in\text{Item}}\sum_{i^{\prime}\in\text{Item}}\sum_{j\in\text{Site}}\sum_{j^{\prime}\in\text{Site}}f_{ii^{\prime}}d_{jj^{\prime}}q_{ij}q_{i^{\prime}j^{\prime}},\end{split}$
(12)

$F_{i}^{\text{Item}}(\mathbf{q})=\sum_{j\in\text{Site}}q_{ij},\quad\forall
i\in\text{Item},$ (13)

$F_{j}^{\text{Site}}(\mathbf{q})=\sum_{i\in\text{Item}}q_{ij},\quad\forall
j\in\text{Site},$ (14)

we can transform the QUBO representation (5) to the following, like (7):

$\begin{array}[]{ll}\text{minimize}&f_{0}(\mathbf{q})+\frac{1}{2}\sum_{i}\lambda\left(F_{i}^{\text{Item}}(\mathbf{q})-1\right)^{2}\\\
&+\frac{1}{2}\sum_{j}\lambda\left(F_{j}^{\text{Site}}(\mathbf{q})-1\right)^{2}\\\
\text{ subject to }&\mathbf{q}\in\\{0,1\\}^{n^{2}}.\end{array}$ (15)

For convenience, we define
$(F_{1},\ldots,F_{2n})=(F_{1}^{\text{Item}},F_{2}^{\text{Item}},\ldots,F_{n}^{\text{Item}},F_{1}^{\text{Site}},F_{2}^{\text{Site}},\ldots,F_{n}^{\text{Site}})$
in (15), obtaining

$\begin{array}[]{ll}\text{minimize}&f_{0}(\mathbf{q})+\frac{1}{2}\sum_{k=1}^{2n}\lambda\left(F_{k}(\mathbf{q})-1\right)^{2}\\\
\text{ subject to }&\mathbf{q}\in\\{0,1\\}^{n^{2}}.\end{array}$ (16)

Similar to the transformation from (7) to (8), we can obtain

$f_{\text{newec}}(\mathbf{q};\mathbf{v})=f_{0}(\mathbf{q})-\sum_{k=1}^{2n}v_{k}F_{k}(\mathbf{q}),$
(17)

using the Hubbard–Stratonovich transformation [9][10], which eliminates the
quadratic term caused by the constraints.

As explained in Section II-C, solutions from the quantum annealer follow a
Boltzmann distribution. In particular, when $f_{\text{newec}}$ is an energy
function, the Boltzmann distribution is

$Q(\mathbf{q})^{t}=Q\left(\mathbf{q};\mathbf{v}^{t}\right)=\frac{\exp\left(-\beta
f_{\text{newec}}\left(\mathbf{q};\mathbf{v}^{t}\right)\right)}{Z(\mathbf{v}^{t})},$
(18)

and its partition function is

$Z(\mathbf{v})^{t}=\sum_{\mathbf{q}}\exp{\left(-\beta
f_{\text{newec}}\left(\mathbf{q};\mathbf{v}^{t}\right)\right)}.$ (19)

To avoid violating the constraints, we update $\mathbf{v}$ as

$v_{k}^{t+1}\leftarrow v_{k}^{t}+\eta^{t}\left(1-\left\langle
F_{k}(\mathbf{q})\right\rangle_{Q^{t}}\right),\quad\forall
k\in\\{1,2,\ldots,2n\\}.$ (20)

To compute $\left\langle F_{k}(\mathbf{q})\right\rangle_{Q^{t}}$, we use the
average $\frac{1}{N}\sum_{i=1}^{N}F_{k}(\mathbf{q}_{i})$ by applying $F$ to
solutions $(\mathbf{q}_{1},\mathbf{q}_{2},\ldots,\mathbf{q}_{N})$ sampled from
the quantum annealer.

In the above, we resolved the difficulty of reducing the problem size by
applying the Ohzeki method to optimization problem (4), which was proposed in
[4]. However, while solutions obtained by the Ohzeki method are all infeasible
in our experiment in Section IV, they satisfy some of the one-hot constraints.
Thus, we can efficiently transform infeasible solutions to feasible ones using
the BFHA, which we propose in the next section.

###  III-B Bit-flip Heuristic in Postprocessing

The Ohzeki method allows solving larger problems with one-hot constraints.
However, when we optimize (4) using the Ohzeki method, we obtain no feasible
solutions, as described in Section IV. Figure 2 shows the number of one-hot
constraints violated in each of 10000 solutions generated by sampling 1000
solutions 10 times with D-Wave Advantage using the Ohzeki method. For example,
a size-8 problem has 16 one-hot constraints. The number of violated one-hot
constraints equals zero in a feasible solution. From Figure 2, there are no
solutions with 0 violations for any problem size, meaning all solutions are
infeasible. Because it is difficult to obtain solutions satisfying all
constraints using the Ohzeki method, we propose the BFHA, which transfers
infeasible solutions violating one-hot constraints, such as (4), into nearby
feasible solutions in terms of the Hamming distance. Before describing the
BFHA, we show that the problem of transferring an infeasible solution to a
nearby feasible solution is the assignment problem where $\bar{q}_{ij}$
represents the infeasible solutions obtained by the Ohzeki method, and
$q_{ij}$ represents the solution after moving by the minimum Hamming distance.
The objective function consists of the sum of $1-q_{ij}$ for $(i,j)$
satisfying $\bar{q}_{ij}=1$ and the sum of $q_{ij}$ for $(i,j)$ satisfying
$\bar{q}_{ij}=0$. In other words, the objective function represents the
Hamming distance between $\bar{q}_{ij}$ and $q_{ij}$. (21) is thus the problem
of finding a feasible solution $q_{ij}$ satisfying the minimum Hamming
distance to $\bar{q}_{ij}$. The Hungarian method [11][12] can similarly solve
this problem as an assignment problem. By solving this optimization problem,
an infeasible solution with constraint violations can be transferred to a
feasible solution with the minimum Hamming distance. However, quantum
annealing usually samples many solutions at once, so a faster algorithm is
needed. We thus propose the BFHA, which is faster than these algorithms.

Figure 2: The numbers of one-hot constraints violated in 10000 solutions
obtained using the Ohzeki method with D-Wave Advantage. Vertical axes show the
number of solutions, and horizontal axes show numbers of violated one-hot
constraints. No problems had feasible solutions for all sizes.

$\begin{array}[]{ll}\text{ minimize
}&\sum_{\bar{q}_{ij}=1}(1-q_{ij})+\sum_{\bar{q}_{ij}=0}q_{ij}\\\ \text{
subject to }&\sum_{i\in\text{Item}}q_{ij}=1,\quad\forall j\in\text{Site},\\\
&\sum_{j\in\text{Site}}q_{ij}=1,\quad\forall i\in\text{Item},\\\
&q_{ij}\in\\{0,1\\},\quad\forall i\in\text{Item},\quad\forall
j\in\text{Site},\end{array}$ (21)

The following briefly describes the BFHA. In the case of problems with
multiple one-hot constraints, such as the QAP and TSP, the solutions can be
arranged as a square matrix. The sum of each row and column is $1$ in a
feasible solution. Also, a solution is feasible if and only if
$\sum_{i=1}^{n}q_{il}+\sum_{j=1}^{n}q_{kj}-2=0,\;\forall k,l$ holds.
$\sum_{i=1}^{n}q_{il}+\sum_{j=1}^{n}q_{kj}-2$ represents difference from the
feasible solution, which is expressed in line 4 of the BFHA (Algorithm 1).
When the input solution to the BFHA is infeasible, the algorithm attempts to
eliminate this gap.

Focusing on lines 6–9 of the algorithm, if $Violation$ has elements with a
value of $1$ or more, the algorithm chooses the $(i,j)$ with the largest
$Violation_{ij}$ satisfying $q_{ij}=1$. Then, the algorithm transforms
$q_{ij}=1$ to $q_{ij}=0$ at $(i,j)$ and updates $Violation$. If no element of
$Violation$ has a value of $1$ or more, and there are elements whose value is
less than or equal to $-1$, then the algorithm executes lines 10-13. The
algorithm chooses the $(i,j)$ with the smallest $Violation_{ij}$ satisfying
$q_{ij}=0$ and transforms $q_{ij}=0$ to $q_{ij}=1$ at $(i,j)$. By repeating
this procedure, an infeasible solution can be shifted to a nearby feasible
solution in terms of the Hamming distance.

Figure 3: Bit-Flip Heuristic Algorithm (BFHA) Figure 4: QAP optimization

To show that the BFHA can always map infeasible solutions to feasible ones, we
first define $X:=[\mathbf{q}_{ij}]$, an $n\times n$ matrix,
$V:=[\mathbf{v}_{ij}]$ where
$\mathbf{v}_{ij}=\sum_{l=1}^{n}q_{il}+\sum_{k=1}^{n}q_{kj}-2$, $U:=X\circ V$
with $\circ$ denoting the Hadamard product, and $D:=(\mathbf{1}-X)\circ V$
where $\mathbf{1}=[1]_{n\times n}$.

###### Theorem 1

The BFHA can always map infeasible solutions to feasible solutions.

###### Proof:

$\mathbf{q}$ is a feasible solution if and only if $V=\mathbf{0}$, and
$V=\mathbf{0}$ if and only if $(U=\mathbf{0})\land(D=\mathbf{0})$. Next, we
show that $\sum_{ij}U_{ij}$ is strictly reduced at lines 6-9 of the algorithm.
We write $U$ as $U^{t}$ to clearly represent the $t$-th loop. If the largest
element of $U^{t}$ is $(i^{\prime},j^{\prime})$, lines 6-9 of the algorithm
are equivalent to setting $U_{i^{\prime}j^{\prime}}^{t}=0$. Thus, using the
fact that $U^{t+1}_{i^{\prime}j^{\prime}}=0\times
v_{i^{\prime}j^{\prime}}^{t+1}=0$, we have
$U^{t}_{i^{\prime}j^{\prime}}-U^{t+1}_{i^{\prime}j^{\prime}}=U^{t}_{i^{\prime}j^{\prime}}\geq
1$. From the above, lines 6-8 of the BFHA strictly reduce $\sum_{ij}U_{ij}$.

We next show that lines 10-13 of this algorithm do not generate elements of
$U$ that are $1$ or more. In other words, we prove that lines 10-13 do not
violate $U=0$. When the smallest element in $D$ satisfies $D_{ij}=-2$, no
element equals $1$ in the $i$-th row and $j$-th column of $X$. Therefore, the
operation at line 12 on $(i,j)$ satisfying $D_{ij}=-2$ has no effect on $U$.We
show that the minimum value of the element of $D$ does not become $-1$ when
$U=0$. We prove this by contradiction. Suppose the value of the smallest
element of $D$ is $-1$. Then for the minimum element $(i,j)$ of $D$, $1$
always exists in the row or column containing $(i,j)$ of $X$. That is, there
are $n$ or more elements in $X$ that equal $1$, but $U=0$ only when the
algorithm is finished, which is a contradiction. From the above, lines 10-13
of this algorithm do not violate $U=0$. This algorithm strictly reduces
$\sum_{ij}U_{ij}$ by lines 6-9 of the algorithm. Thus, $U=0$ will always be
true. For the same reason, $D=0$ is always satisfied by lines 10-13. Also,
$U=0$ is not violated by lines 10-13 of this algorithm, thus establishing
$(U=\mathbf{0})\land(D=\mathbf{0})$. Therefore, the BFHA can always map
infeasible solutions to feasible ones. ∎

As mentioned above, it is difficult to obtain feasible solutions using only
the Ohzeki method. We thus propose Algorithm 2, which always yields feasible
solutions by applying both the Ohzeki method and the BFHA. Specifically,
parameter $v$ is updated according to the Ohzeki method, and the BFHA is
applied to the obtained solution every time, preserving the solution with the
lowest objective function value. We can thereby solve problems (4) with sizes
5 to 20 using D-Wave Advantage, as described in Section IV.

##  IV RESULT AND DISCUSSION

###  IV-A BFHA Performance

We conducted experiments to demonstrate the BFHA speed and its ability to
transfer infeasible solutions to nearby feasible solutions. The computing
environment was an Intel Core i7-7700HQ 2.80 GHz CPU with 4 cores and 8
threads. Using the Ohzeki method, we generate 1000 infeasible solutions for
each size of problem (4). We measured the following two indicators:

  * •

Computation time (in seconds) per 1000 infeasible solutions

  * •

Hamming distance between an infeasible solution and a transferred feasible
solution

For the above two indices, we conducted 10 experiments using the following
three methods:

  * •

Solve (21) using Gurobi [13]

  * •

Solve (21) by the Hungarian method [11][12], using the Munkres implementation
for Python [14]

  * •

Use the BFHA

Table I summarizes mean computation times and errors in Hamming distance, with
“Gurobi time [$\mathrm{s}$],” “Hungarian time [$\mathrm{s}$],” and “Bitflip
time [$\mathrm{s}$]” respectively showing computation times in seconds under
the Gurobi, Hungarian, and BFHA methods. The “Gap [$\mathrm{\char 37\relax}$]”
column shows error between the minimum Hamming distance and that under the
BFHA, indicating high accuracy. This table also confirms that the BFHA is much
faster than the other algorithms, and that the increase in computation time
for each size is larger with the Hungarian method than with the BFHA. However,
the Gurobi computation time appears to increase at a slower rate, because it
appears to take a long time to formulate the problem, thus gradually
increasing calculation times. After building Violation, the BFHA performs only
simple calculations such as adding and subtracting elements in matrix rows and
columns. The BFHA can reduce the number of iterations and can find a feasible
solution more quickly when there are few one-hot constraints violated among
the infeasible solutions. Since some one-hot constraints are satisfied using
the Ohzeki method, the BFHA is considered to be fast at finding a feasible
solution.

TABLE I: Comparison of speed and Hamming distance Size | Gurobi time [$\mathrm{s}$] | Hungarian time [$\mathrm{s}$] | Bitflip time [$\mathrm{s}$] | Gap [$\mathrm{\char 37\relax}$]
---|---|---|---|---
5 | $3.991\times 10^{0}$ | $9.516\times 10^{-2}$ | $3.973\times 10^{-3}$ | $1.437$
6 | $4.953\times 10^{0}$ | $1.605\times 10^{-1}$ | $4.476\times 10^{-3}$ | $2.256$
7 | $5.560\times 10^{0}$ | $1.962\times 10^{-1}$ | $4.963\times 10^{-3}$ | $1.371$
8 | $6.122\times 10^{0}$ | $2.236\times 10^{-1}$ | $5.886\times 10^{-3}$ | $1.712$
9 | $6.743\times 10^{0}$ | $2.913\times 10^{-1}$ | $7.619\times 10^{-3}$ | $2.874$
10 | $7.385\times 10^{0}$ | $3.393\times 10^{-1}$ | $7.836\times 10^{-3}$ | $2.588$
11 | $8.096\times 10^{0}$ | $4.695\times 10^{-1}$ | $1.057\times 10^{-2}$ | $3.816$
12 | $9.010\times 10^{0}$ | $5.171\times 10^{-1}$ | $1.167\times 10^{-2}$ | $4.159$
13 | $9.784\times 10^{0}$ | $6.279\times 10^{-1}$ | $1.323\times 10^{-2}$ | $4.517$
14 | $1.061\times 10^{1}$ | $6.930\times 10^{-1}$ | $1.404\times 10^{-2}$ | $4.596$
15 | $1.173\times 10^{1}$ | $8.518\times 10^{-1}$ | $1.790\times 10^{-2}$ | $6.233$
16 | $1.242\times 10^{1}$ | $1.061\times 10^{0}$ | $2.161\times 10^{-2}$ | $6.844$
17 | $1.162\times 10^{1}$ | $1.117\times 10^{0}$ | $2.881\times 10^{-2}$ | $7.491$
18 | $1.375\times 10^{1}$ | $1.151\times 10^{0}$ | $3.570\times 10^{-2}$ | $7.395$
19 | $1.331\times 10^{1}$ | $1.237\times 10^{0}$ | $4.421\times 10^{-2}$ | $6.712$
20 | $1.454\times 10^{1}$ | $1.434\times 10^{0}$ | $4.702\times 10^{-2}$ | $6.924$

###  IV-B Larger QAP Optimization on an E-commerce Website

####  IV-B1 Problem setting

We randomly created 10 problems with 5 to 20 items and formulated them as (4).
We set $w=0.5$ in (4), and set $\eta^{t}=0.1\times 0.92^{t-1}$ in Algorithm 2
(except that $\eta^{t}=\eta^{9}$ if $t\geq 9$). We set the maximum number of
iterations of Algorithm 2 to 30. We solved each problem one at a time using
Algorithm 2. We used Advantage_system1.1 with 5436 qubits, which is the latest
version of D-Wave Advantage as of November 2020. We set quantum annealer
parameters as num_reads = 1000, annealing time = 50
$\mathrm{\SIUnitSymbolMicro s}$. We used the method of [15] to embed the
problem in the quantum annealer.

####  IV-B2 Result

Each problem was solved once by three methods: Algorithm 2, only the Ohzeki
method, and a method that directly solves (5). Table II summarizes means of
objective function values for the 10 best solutions by each method. Note that
the Ohzeki method alone obtained no feasible solutions after 30 times with
1000 samples, so Table II does not list them. The “With bitflip” and “QUBO”
columns respectively show the results of Algorithm 2 and the direct method.
“Score” and “Time [$\mathrm{s}$]” are respectively means of the objective
function value for the best obtained solution and the computation time, and
“Opt” is the average of objective function values of the exact solutions.
Cases where no feasible solution was obtained are denoted as “infeasible,” and
problems of a size that cannot be embedded in the quantum annealer graph are
denoted as “cannot solve.”

We define the computation time for Algorithm 2 as the sum of the time for
annealing, the BFHA execution time, the time for calculating expectations, and
the time for updating parameters. Communication times between Japan and D-Wave
Systems’ quantum annealer in Canada is not included in the computation time.
The computation time for “QUBO” includes only the time for quantum annealing.
Calculations are terminated if an exact solution is found. Table II shows that
Algorithm 2 (“With bitflip”) obtains overall good results. However, the
difference between the obtained scores and “Opt” increases with size. For size
20 in particular, the solution accuracy deteriorated immediately. We found no
feasible solutions by the Ohzeki method alone. We found no feasible solution
for 6 out of 10 size-7 problems, so the average objective function value of
the remaining four best solutions is given as the “QUBO” value. Similarly, we
found no feasible solution for 8 out of 10 size-8 problems, and so use the
remaining two solutions. We could obtain no feasible solutions for size 9 or
more. From this, we can see that it is difficult to solve a problem with
multiple one-hot constraints, such as the QAP, in the QUBO representation (5).
In terms of speed, the computation times of Algorithm 2 show a gradual
increase. The computation time of the quantum annealer does not change, and
the increase in the BFHA computation time is not large. Therefore, the
increase in calculation time is mostly due to line 8 in Algorithm 2.

TABLE II: Results of a numerical experiment with D-Wave Advantage.  | With bitflip | QUBO |
---|---|---|---
Size | Score | Time [$\mathrm{s}$] | Score | Time [$\mathrm{s}$] | Opt
5 | $-2.884\times 10^{0}$ | $1.165\times 10^{-1}$ | $-2.870\times 10^{0}$ | $7.200\times 10^{-1}$ | $-2.884\times 10^{0}$
6 | $-3.838\times 10^{0}$ | $2.856\times 10^{-1}$ | $-3.588\times 10^{0}$ | $1.385\times 10^{0}$ | $-3.838\times 10^{0}$
7 | $-4.779\times 10^{0}$ | $9.658\times 10^{-1}$ | $-2.599\times 10^{0}$ | $1.500\times 10^{0}$ | $-4.812\times 10^{0}$
8 | $-5.351\times 10^{0}$ | $2.901\times 10^{0}$ | $-1.944\times 10^{0}$ | $1.500\times 10^{0}$ | $-5.530\times 10^{0}$
9 | $-5.949\times 10^{0}$ | $4.122\times 10^{0}$ | infeasible | infeasible | $-6.245\times 10^{0}$
10 | $-6.308\times 10^{0}$ | $4.003\times 10^{0}$ | infeasible | infeasible | $-6.906\times 10^{0}$
11 | $-6.964\times 10^{0}$ | $5.904\times 10^{0}$ | infeasible | infeasible | $-7.880\times 10^{0}$
12 | $-7.307\times 10^{0}$ | $4.193\times 10^{0}$ | infeasible | infeasible | $-8.612\times 10^{0}$
13 | $-7.876\times 10^{0}$ | $4.566\times 10^{0}$ | cannot solve | cannot solve | $-9.629\times 10^{0}$
14 | $-8.231\times 10^{0}$ | $4.934\times 10^{0}$ | cannot solve | cannot solve | $-1.040\times 10^{1}$
15 | $-8.778\times 10^{0}$ | $5.261\times 10^{0}$ | cannot solve | cannot solve | $-1.133\times 10^{1}$
16 | $-9.169\times 10^{0}$ | $6.187\times 10^{0}$ | cannot solve | cannot solve | $-1.211\times 10^{1}$
17 | $-9.482\times 10^{0}$ | $6.593\times 10^{0}$ | cannot solve | cannot solve | $-1.307\times 10^{1}$
18 | $-9.915\times 10^{0}$ | $7.237\times 10^{0}$ | cannot solve | cannot solve | $-1.156\times 10^{1}$
19 | $-1.042\times 10^{1}$ | $7.853\times 10^{0}$ | cannot solve | cannot solve | $-1.209\times 10^{1}$
20 | $-1.078\times 10^{1}$ | $9.857\times 10^{0}$ | cannot solve | cannot solve | $-1.563\times 10^{1}$

####  IV-B3 Discussion

The proposed method, Algorithm 2, showed good results overall. There are two
possible reasons for this. First, by removing the one-hot constraint,
connections between logical variables become sparse, so we can transform the
original problem into one that is easy for the quantum annealer to solve. In
[16], D-Wave 2000Q was used to solve the max-cut problem for regular graphs
with different degrees. As the degree increases, that is, as the coupling
becomes denser, the problem becomes more difficult to solve. Because $d$ in
(4) makes problems sparse and we remove the one-hot constraints, we can
transform the original problem into an easier problem for the quantum
annealer. Secondly, the BFHA transfers infeasible solutions produced by the
quantum annealer to nearby feasible ones according to Hamming distance.

However, the accuracy of size-20 problems was worse, and the same trend is
observed in size-21 and above problems. We assume there are three main reasons
for this. First, the chain is long and the number of qubits used is increased.
We find that as the problem size increases, the chain length and the number of
used qubits increase in the Pegasus graph. A chain is qubits used to represent
a logical variable. Table III shows the average number of couplers and the
average chain length for each problem size in the D-Wave Advantage. We can see
that as the problem size increases, the chain length and number of required
qubits increases. As reported in [16], solution accuracy is known to becomes
worse as the number of network qubits increases, which may have affected the
results. Secondly, quantum annealer noise may have prevented correct sampling.
It is known that more quantum annealer noise worsens solution accuracy [17],
and this occurs with D-Wave Systems’ quantum annealer. The third possibility
is that we may not be able to obtain a good solution with the Ohzeki method,
which is a heuristic algorithm with no theoretical guarantees, so the
algorithm may have deteriorated solution accuracy. Further analysis of this
point is a subject for future study.

###  IV-C Contribution

Our study made two main contributions. First, we proposed a heuristic
algorithm named the bit-flip heuristic algorithm (BFHA), which moved solutions
violating the one-hot constraint to nearby feasible solutions in terms of the
Hamming distance. We also showed that the algorithm always transferred
infeasible solutions to feasible ones. Furthermore, we numerically confirmed
that the algorithm was fast and had high accuracy, meaning a small Hamming
distance between the feasible solution assigned by the algorithm and the
original infeasible solution. Second, we proposed a method for solving the QAP
with quantum annealing by applying the BFHA to solutions obtained by the
Ohzeki method. In a numerical experiment, we applied the proposed method to
solve a sparse QAP that was used in areas such as item listing on an
E-commerce website [4]. Nishimura et al. [4] solved this problem with up to 8
items using D-Wave 2000Q. If D-Wave Advantage had been used, problems with up
to 12 items could have been solved. In this study, we successfully solved a
QAP of size 19 with high accuracy for the first time using a quantum annealer
without splitting.

TABLE III: Numbers of logical variables, required qubits, and required couplers, and average chain length for each problem size in the D-Wave Advantage. Size | Logical variables | Required Qubits | Required Couplers | Average chain length
---|---|---|---|---
5 | 25 | 32 | 80 | 1.280
6 | 36 | 53 | 150 | 1.472
7 | 49 | 83 | 252 | 1.694
8 | 64 | 123 | 392 | 1.922
9 | 81 | 176 | 576 | 2.173
10 | 100 | 248 | 810 | 2.480
11 | 121 | 325 | 1100 | 2.686
12 | 144 | 412 | 1452 | 2.861
13 | 169 | 529 | 1872 | 3.130
14 | 196 | 634 | 2366 | 3.235
15 | 225 | 784 | 2940 | 3.484
16 | 256 | 958 | 3600 | 3.742
17 | 289 | 1223 | 4352 | 4.232
18 | 324 | 1472 | 5202 | 4.543
19 | 361 | 1710 | 6156 | 4.737
20 | 400 | 2001 | 7220 | 5.003

##  V CONCULUSION

In this study, we obtained approximate solutions for a sparse QAP of size 19
using D-Wave Advantage. We also proposed the BFHA, a method for moving
infeasible solutions to a nearby feasible solution in terms of the Hamming
distance. We also proved that the algorithm always moves infeasible solutions
to feasible solutions, and numerically confirmed that the algorithm can move
to the nearest feasible solution. We confirmed that the algorithm is faster
than the compared algorithms. In a practical use case, we optimized item
placement on an E-commerce website because this sparse QAP has been used in
areas such as item listing on an E-commerce website. In [4], D-Wave 2000Q
solved problems of up to size 8. When using the same method with D-Wave
Advantage, it is impossible to solve problems of size 13 or larger. In
addition, the Ohzeki method proposed in [7] provided no feasible solutions to
our problem. In contrast, using the BFHA with the Ohzeki method allowed
optimization of up to 19 items in D-Wave Advantage with good accuracy.

Our proposed method is applicable to the real world applications such as item
listing on an E-commerce website[4] and minimizing the backtrack of jobs in
the Generalized Flowline (GFL)[18], where GFL is a flow line on which jobs
flow downstream and the backtrack is the movement of the auto guided vehicle
(AGV) against the ideal direction of GFL. Moreover, many real-world
optimization problems involve one-hot constraints, so the proposed method has
a wide range of applications, not only to the QAP but also to other problems
with one-hot constraints. In this sense, the proposed method can thus make
significant contributions to the development of real-world applications using
quantum annealers. We obtained approximate solutions to size-19 QAP by the
proposed method, but it is known that such approximate solutions can be
quickly obtained even on classical computers. Like other studies [3][4], this
method is therefore not practical at present. Even so, quantum annealer
technology is rapidly advancing, from 128 qubits in 2011 to 5640 qubits in
2020. In addition, noise has been reduced and more accurate solutions can be
found, and quantum annealer hardware will continue to evolve in the future.
For this reason, it is possible that this method will be practical for real-
world applications in the future. However, it is necessary in future studies
to address the problem of worsening accuracy for problems of size 20 and
above.

## References

  * [1] T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse ising model,” _Physical Review E_ , vol. 58, no. 5, pp. 5355–5363, Nov. 1998, DOI: 10.1103/PhysRevE.58.5355.
  * [2] D-Wave Systems. 2019. “Next-generation topology of d-wave quantum processors.” D-Wave White Papers. https://www.dwavesys.com/sites/default/files/14-1026A-C_Next-Generation-Topology-of-DW-Quantum-Processors.pdf. (accessed Dec. 3, 2020).
  * [3] F. Neukart, G. Compostella, C. Seidel, D. V. Dollen, S. Yarkoni, and B. Parney, “Traffic flow optimization using a quantum annealer,” _Frontiers in ICT_ , vol. 4, pp. 29, Dec. 2017, DOI: 10.3389/fict.2017.00029.
  * [4] N. Nishimura, K. Tanahashi, K. Suganuma, M. J. Miyama, and M. Ohzeki, “Item listing optimization for e-commerce websites based on diversity,” _Frontiers in Computer Science_ , vol. 1, pp. 2, July 2019, DOI: 10.3389/fcomp.2019.00002.
  * [5] T. C. Koopmans and M. Beckmann, “Assignment problems and the location of economic activities,” _Econometrica_ , vol. 25, no. 1, pp. 53–76. Jan. 1957, DOI: 10.2307/1907742.
  * [6] D-Wave Systems. 2017. “Partitioning optimization problems for hybrid classical/quantum execution.” D-Wave White Papers. https://docs.ocean.dwavesys.com/projects/qbsolv/en/latest/_downloads/bd15a2d8f32e587e9e5997ce9d5512cc/qbsolv_techReport.pdf. (accessed Dec. 3, 2020).
  * [7] M. Ohzeki, “Breaking limitation of quantum annealer in solving optimization problems under constraint,” _Scientific Reports (Nature Publisher Group)_ , vol. 10, no. 1, Feb. 2020.
  * [8] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” _science_ , vol. 220, no. 4598, pp. 671–680, DOI: 10.1126/science.220.4598.671.
  * [9] J. Hubbard, “Calculation of partition functions,” _Physical Review Letters_ , vol. 3, no. 2, pp. 77–78, July. 1959, DOI: 10.1103/PhysRevLett.3.77.
  * [10] R. L. Stratonovich, “On a method of calculating quantum distribution functions,” _SPhD_ , vol. 2, pp. 416, 1957.
  * [11] H. W. Kuhn, “The hungarian method for the assignment problem,” _Naval research logistics quarterly_ , vol. 2, no. 1–2, pp. 83–97, 1955, DOI: 10.1002/nav.3800020109.
  * [12] J. Munkres, “Algorithms for the assignment and transportation problems,” _Journal of the society for industrial and applied mathematics_ , vol. 5, no. 1, pp. 32–38, 1957, DOI: 10.1137/0105003.
  * [13] _gurobi_. (Gurobi 9.1), Gurobi Optimization. Accessed: Dec. 3, 2020. [Online]. Available: http://www.gurobi.com
  * [14] _munkres — Munkres implementation for Python_. (1.1.2), B. M. Clapper. Accessed: Dec. 3, 2020. [Online]. Available: http://software.clapper.org/munkres/index.html
  * [15] J. Cai, W. G. Macready, and A. Roy, “A practical heuristic for finding graph minors,” 2014. [Online]. Available: arXiv:1406.2741 [quant-ph].
  * [16] R. Hamerly _et al._ , “Experimental investigation of performance differences between coherent ising machines and a quantum annealer,” _Science advances_ , vol. 5, no. 5, eaau0823, May. 2019, DOI: 10.1126/sciadv.aau0823.
  * [17] D-Wave Systems. 2019. “Improved coherence leads to gains in quantum annealing performance,” D-Wave White Papers. https://www.dwavesys.com/sites/default/files/14-1037A-A_Improved_coherence_leads_to_gains_QA_performance.pdf (accessed Dec. 3, 2020).
  * [18] D. Gong, G. Yamazaki, M. Gen, and W. Xu, “A genetic algorithm method for one-dimensional machine location problems,” _International Journal of Production Economics_ , vol. 60, pp. 337–342, 1999, DOI: 10.1016/S0925-5273(98)00163-7.

000All company or product names mentioned herein are trademarks or registered
trademarks of their respective owners.