# Diversity metric for evaluation of quantum annealing

Alex Zucca  D-Wave Systems Inc.  Simon Fraser University  Hossein Sadeghi
D-Wave Systems Inc.  Masoud Mohseni  Google Inc.  Mohammad H. Amin  D-Wave
Systems Inc.  Simon Fraser University

###### Abstract

Solving discrete NP-hard problems is an important part of scientific
discoveries and operations research as well as many commercial applications. A
commonly used metric to compare meta-heuristic solvers is the time required to
obtain an optimal solution, known as time to solution. However, for some
applications it is desirable to have a set of high-quality and diverse
solutions, instead of a single optimal one. For these applications, time to
solution may not be informative of the performance of a solver, and another
metric would be necessary. In particular, it is not known how well quantum
solvers sample the configuration space in comparison to their classical
counterparts. Here, we apply a recently introduced collective distance measure
in solution space to quantify diversity by Mohseni, et. al. Mohseni et al.
(2021) and, based on that, we employ time-to-diversity as a metric for
evaluation of meta-heuristics solvers. We use this measure to compare the
performance of the D-Wave quantum annealing processor with a few classical
heuristic solvers on a set of synthetic problems and show that D-Wave quantum
annealing processor is indeed a competitive heuristic, and on many instances
outperforms state-of-the-art classical solvers, while it remains on par on
other instances. This suggests that a portfolio solver that combines quantum
and classical solutions may win over all solvers.

##  I Introduction

The development of quantum annealing (QA) processors (Johnson et al., 2011),
has sparked a renewed interest on combinatorial optimization for industrial
problems. The potential advantage of QA processors lies in the fact that they
have access to quantum resources, such as superposition and entanglement among
variables, which could lead to the possibility of non-trivial quantum
tunneling or collective moves in configuration space that are absent in local
stochastic classical heuristics. The quantum annealer is usually tested
against some classical algorithms, such as simulated annealing (SA), parallel
tempering (PT), and parallel tempering with iso-energetic cluster moves
(PT+ICM) Zhu et al. (2015, 2016), with the general trend that QA competes
fairly well with the state of the art heuristics. Due to the limitations of
the current QA processors in terms of number of variables and their inter-
connectivity, the search for quantum speed-up is usually performed on a set of
synthetic benchmark problems Hen et al. (2015); King et al. (2015a, 2017);
Albash and Lidar (2018); Mandrà and Katzgraber (2018), where the “hardness” of
the problems is, in some cases, tuneable via one or more parameters. However,
despite many attempts to detect quantum speedup in the current generation of
quantum annealers over classical algorithms Dickson et al. (2013); Boixo et
al. (2014); Katzgraber et al. (2015); Heim et al. (2015); Katzgraber et al.
(2014, 2015); Heim et al. (2015); Hen et al. (2015); Albash et al. (2015);
Martin-Mayor and Hen (2015); King et al. (2015a); Marshall et al. (2016);
Denchev et al. (2016); King et al. (2017); Albash and Lidar (2018), the
existence of an improved scaling in optimization still remains elusive,
although such a scaling advantage has recently been reported in quantum
simulation of magnetic materials King et al. (2021).

Most previous studies reported the performances in terms of time-to-solution
(TTS), typically defined as the total time (which could be the sum of several
repetitions) required to reach an optimal solution with 99% certainty. It is
clear that in order to calculate TTS, knowledge of the optimal solution is
necessary, which can become intractable for large problems. Moreover, TTS is
defined based on obtaining a single solution and therefore does not give any
information on how widely the returned solutions are distributed in the
configuration space. Other metrics such as the Time-To-Target (TTT) introduced
in King et al. (2015b) are designed to provide a benchmarking quantity that
does not become intractable as the number of variables increases. This is
achieved by considering sub-optimal solutions provided by a reference solver
(e.g. the quantum annealer) and estimating how long other solvers take to
achieve the same level of sub-optimality. Other studies instead, analyzed the
sampling abilities of QA. In particular it was noted in Matsuda et al. (2009)
that transverse field driver Hamiltonians like the one implemented in current
quantum annealers result in an unfair sampling of degenerate ground states.
Indeed this was confirmed in a systematic study of the D-Wave 2X annealer
Mandrà et al. (2017) and it was suggested the more complex driver Hamiltonian
may help reduce the bias, although with the conditions of being dense enough
Könz et al. (2018). Nevertheless, QA can be complementary to classical
algorithms to obtain “classically rare” solutions Zhang et al. (2017).
Recently, fair sampling performances of state-of-the-art classical algorithms
such as PT and PT+ICM were put under scrutiny Zhu et al. (2019). Although it
emerged that PT+ICM performs better than PT, it was also noted the
performances degrade as the Hamming distance between ground states increase. A
different point of view was taken in King et al. (2019) where the authors
defined the Time-To-All-Valleys (TTAV) metric to measure the ability of a
solver to obtain at least one sample in an island of ground states, where the
latter is defined as a configuration of ground states that are connected in
the logical Hamming space. Since moving within the same valley can be done
with minimal post-processing, being able to reach quickly all the valleys is
of critical importance.

In this work we adopt a similar point of view, but we do not restrict to
diverse ground state configurations. Rather, we focus on the ability of
solvers to reach diverse sub-optimal solutions. In fact, due to noise,
calibration errors, and finite temperature, QA processors are expected to
produce a pool of diverse, high-quality but sub-optimal solutions. Diversity
of high-performing population is of crucial importance in genetic programming,
evolutionary algorithms, quality-diversity optimization Pugh et al. (2016),
contingency planning, robust optimization, and machine learning applications
Ulrich et al. (2010); Ulrich and Thiele (2011); Gao et al. (2015); Alexander
et al. (2017); Neumann et al. (2018). Recently, a new diversity measure was
introduced based on the number of high-quality solutions, within a target
approximation ratio, that have large mutual Hamming distance larger than a
certain threshold Mohseni et al. (2021). Here, we adopt such diversity measure
to quantify the performance of the D-Wave 2000Q annealer in terms of time-to-
diversity (TTD) - the time the solver takes to return a sample set with target
diversity.

The paper is structured as follows. In Sect. II we review how the D-Wave
quantum processing unit (QPU) solves an Ising minimization problem. In Sect.
III we review the Hamming-distance-based _diversity measure_ that can be used
to characterize the ability of a solver to return diverse samples Mohseni et
al. (2021). Similarly to the definition of TTS, we employ the concept of TTD
and we introduce two ways of obtaining its upper and lower bounds. We present
and discuss the results in Sect. IV. We then conclude in Sect. V.

##  II Overview of the D-Wave hardware

The D-Wave quantum annealer is designed to solve Ising-like binary
optimization problems, and implements, in the close system approximation, the
following Hamiltonian,

$\displaystyle H$
$\displaystyle=-A(s)\sum_{i}\sigma_{i}^{x}+B(s)\mathcal{H}_{P},$ (1)
$\displaystyle\mathcal{H}_{P}$
$\displaystyle=\sum_{i}h_{i}\sigma_{i}^{z}+\sum_{i,j}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z},$
(2)

where $h_{i}$ and $J_{ij}$ are local and pairwise parameters, respectively,
and $\sigma_{i}^{x}$ ($\sigma_{i}^{z}$) represents the Pauli matrices in $x$
($z$) direction. The set of indices $(i,j)$ forms the adjacency matrix of the
Chimera graph, the architecture of the D-Wave 2000Q quantum annealer Bunyk et
al. (2014). The two functions $A(s)$ and $B(s)$ – known as the annealing
schedule – have a pre-defined time dependent form expressed in terms of the
normalized time $s=t/t_{a}$, with $t_{a}$ being the total anneal time. A
problem of interest is usually mapped (in polynomial time) Lucas (2014) to the
problem of minimizing the energy of an Ising Hamiltonian of the form
above.111While mapping a NP-hard problem to an Ising Hamiltonian is
polynomial, finding the optimal embedding onto the Chimera graph is NP-hard
itself. More time-efficient heuristic ways to find an embedding are used
instead Boothby et al. (2016); Cai et al. (2014). In forward annealing, the
system is initially prepared in the ground state of the Hamiltonian when
$A(0)\gg B(0)$. The energy scale of the driver Hamiltonian – $A(s)$ – is then
slowly reduced, while the energy scale of the problem Hamiltonian – $B(s)$ –
is increased to give, at end of the anneal, $B(1)\gg A(1)\approx 0$. At $s=1$,
the Hamiltonian (1) reduces to the classical energy of the Ising model, at
which point we replace the Pauli matrices with the values of spin represented
by $s_{i}\in\\{\pm 1\\}$,

$\mathcal{H}_{P}=\sum_{i}h_{i}s_{i}+\sum_{i,j}J_{ij}s_{i}s_{j}.$ (3)

If the process is done slowly enough, due to the adiabatic theorem Ehrenfest
(1916); Born and Fock (1928), the system will remain in the ground state of
the instantaneous Hamiltonian, and will return at the end of the anneal the
spins configuration that minimizes the Ising Hamiltonian (3).

In practice, however, the evolution is usually non-adiabatic and is done in
the presence of an environment with finite temperature and Hamiltonian
parameter specification error. As such, instead of obtaining the ground state
configuration at every anneal, one rather obtains a distribution of
configurations with a probability $p$ to obtain a ground state solution. Here
we turn our attention on a benchmarking quantity that does not involve the
achievement of the optimal solution. Rather, we assess the ability to explore
sub-optimal solutions in a different manner. While estimating the solutions to
a sub-optimal point is proven to be polynomial in time Saket (2013), the pre-
factor of such algorithms grows exponentially with the inverse of normalized
residual energy, turning these methods impractical for solving problems of our
interest.

##  III Diversity measure based on mutual Hamming distance

In this Section, we define the _diversity measure_ based on the mutual Hamming
distance between different samples Mohseni et al. (2021), and adopt the
related idea of Time-To-Diversity (TTD). Because the computation of diversity
becomes quickly intractable with the system size, we identify some more
efficient strategies to assess the diversity of the samples returned by
solvers.

Let us consider the set of $2^{N}$ configurations of a binary optimization
problem with Hamiltonian (3) with $N$ being the total number of spins, and let
us denote the ground state energy as $E_{0}$. To restrict our attention to
sub-optimal configurations, for a given solution ${\bf s}^{a}=\\{s^{a}_{i}\\}$
with energy $E({\bf s}^{a})$ we define approximation ratio as

$r({\bf s}^{a})=\frac{E({\bf s}^{a})-E_{0}}{E_{\rm max}-E_{0}},$ (4)

where $E_{\rm max}$ is the maximum energy of the problem Hamiltonian. We
choose a threshold value $0\leq\alpha\ll 1$ for the approximation ratio and a
corresponding energy threshold

$E_{\alpha}=E_{0}+\alpha(E_{\rm max}-E_{0}).$ (5)

Let us define $\mathcal{S}=\\{{\bf s}^{1},\dots,{\bf s}^{M}\\}$ as a set of
$M$ spin configurations ${\bf s}^{a}$ having $E_{0}\leq E({\bf s}^{a})\leq
E_{\alpha}$. For all states ${\bf s}^{a}\in\mathcal{S}$, we can build the
Hamming distance matrix

$d_{ab}=\frac{1}{2}\sum_{i=1}^{N}|s^{a}_{i}-s^{b}_{i}|,$ (6)

where $a$ and $b$ are the indices of solutions and $i$ is the index of
individual binary units (qubits). We then construct a graph
$\mathcal{G}_{\mathcal{S}R}=(\mathcal{V}=\mathcal{S},\mathcal{E})$ where each
solution ${\bf s}^{a}\in\mathcal{S}$ is a node, and the edges are defined only
for the pairs of solutions ${\bf s}^{a},{\bf s}^{b}$ that have Hamming
distance smaller than a normalized radius $R$ through $d_{ab}\leq RN$. An
independent set $V_{\rm IS}$ is a set of vertices, no pair of which are
connected by an edge. The _diversity_ $D_{\mathcal{S}R}$ is defined as the
_independence number_ of the graph $\mathcal{G}_{\mathcal{S}R}$, i.e. the size
of the maximum independent set (MIS)

$D_{\mathcal{S}R}=\max_{V_{\rm IS}\subseteq\mathcal{V}}|V_{\rm IS}|.$ (7)

Note that $E_{\alpha}$ is chosen based on the approximation ratio that
requires knowledge of the ground state energy. However, it is possible to
choose the energy threshold directly (e.g., based on the application) without
knowing the ground state and the following discussions remain unaffected.
Exact calculation of $D_{\mathcal{S}R}$ becomes computationally intractable
when the number of solutions in $\mathcal{S}$ is large. Nevertheless, we can
compute approximated diversity as we discuss the following.

###  III.1 Approximated Diversity Measures

Here, we show that upper and lower bounds on the diversity measure can be
computed with polynomial algorithms. For the lower bound,
$D_{\mathcal{S}R}^{\rm low}$ we use a recursive algorithm that we call _Left
Neighbors Aggregation_ (LNA) algorithm. The algorithm sequentially scans every
vertex and add it to the current independent set if the vertex is not
connected to the previous members of the set. Naturally, the order with which
the samples are scanned is important, and to tighten the bound we analyze many
random orders for the set $\mathcal{S}$. The algorithm is described in details
in App. A.1. The LNA algorithm complexity is $\mathcal{O}(M\cdot
D_{\mathcal{S}R}^{\rm low})$ . In practical cases $D_{\mathcal{S}R}^{\rm
low}\ll M$, so that the LNA scales almost linearly with the number of samples
$M$.

The upper bound $D_{\mathcal{S}R}^{\rm up}$ is computed by approximating the
_chromatic number_ $\chi$ of the complement graph
$\bar{\mathcal{G}}_{\mathcal{S}R}$ 222The chromatic number $\chi(\mathcal{G})$
of a graph $\mathcal{G}$ is the minimum $k$ for which $\mathcal{G}$ admits a
_legal_ coloring of it vertices. A legal coloring of a graph is an assignment
of colors to its vertices so that no pair of adjacent vertices has the same
color. The chromatic number is also an upper bound on the clique number
$\omega(\mathcal{G})$ of the same graph, which is equivalent to the
independence number $\alpha(\bar{G})$ of the complement graph
$\bar{\mathcal{G}}$,
$\alpha(\bar{\mathcal{G}})=\omega(\mathcal{G})\leq\chi(\mathcal{G})$. with a
greedy coloring algorithm Kosowski and Manuszewski (2004); Matula and Beck
(1983). Computing the approximated chromatic number with the greedy algorithm
is $\mathcal{O}(M^{2})$ in the number of samples $M$, as one needs to compute
the distance matrix $d_{ab}$ in advance. The greedy coloring algorithm is
described in App. A.2.

As we show in App. A.3, the lower bound obtained by the LNA algorithm is
tighter than the upper bound and scales similarly as the independence number.
Thus, in the following we will exclusively refer to diversity as its lower
bound computed with the LNA algorithm.

###  III.2 Time-to-Target-Diversity

Inspired by King et al. (2015b) we introduce a strategy to compare the
diversity performance of different solvers. We consider one of the solvers as
a reference solver $\kappa_{r}$ \- in our case the D-Wave QPU - and define a
target diversity in the following way. We obtain a fixed number of samples
from the reference solver, then remove all samples with energy above the
threshold $E_{\alpha}$. This new set of samples, denoted as
$\mathcal{S}_{\alpha}^{r}$, is then used to calculate the diversity. For a
fixed approximation ratio $\alpha$ and hamming distance radius $R$, the target
diversity is defined as $\hat{D}\equiv D_{\mathcal{S}_{\alpha}^{r}R}$.

Consider a set of solutions $\mathcal{S}_{\kappa}(t)$ with energies below
$E_{\alpha}$, returned by the solver $\kappa$ within finite time
$t=t_{a}n_{\rm runs}$, where $n_{\rm runs}$ is the number of repetitions and
$t_{a}$ is the time it takes to run the algorithm once. For simplicity of
notation, we denote the diversity of this set as $D_{\kappa}(t)\equiv
D_{\mathcal{S}_{\kappa}(t)R}$. Let $p_{\kappa}$ be the probability of the
solver $\kappa$ reaching the target diversity, $D_{\kappa}(t)\geq\hat{D}$,
after time $t$. Similarly to TTS, we define Time-to-Target-Diversity (TTD) as
the time it takes for the solver to reach the target diversity $\hat{D}$ with
$99\%$ confidence,

${\rm TTD}_{\kappa}=n_{\rm runs}t_{a}\frac{\log(0.01)}{\log(1-p_{\kappa})}.$
(8)

##  IV Results and discussion

In this Section we compare the DW2000Q quantum annealer against two commonly
used classical solvers, namely Parallel Tempering (PT) and Parallel Tempering
with iso-energetic cluster moves (PT+ICM). In our benchmark we consider the
reference diversity $\hat{D}$ using 50k samples obtained by the DW2000Q
solver.

We begin our analysis using an approximation ratio $\alpha=0.01$ and a radius
factor $R=0.25$. Smaller approximation ratios require significant computation
times. We also expect the performance of the QPU to degrade due to noise and
calibration error. Nevertheless, we studied $\alpha=0.005$ with smaller number
of problems and observed similar behavior as shown in Sect. IV.6.

###  IV.1 Input Classes

We analyzed ${\rm TTD}$ for three different input classes. Each of them is
randomly generated on the underlying Chimera graph of the available D-Wave
QPU. The inputs were generated using only the first $L\times L$ unit cells of
the Chimera graph with $L\in\\{8,12,16\\}$ with 512, 1152 and 2048 qubits
respectively. All the classes have linear terms $h_{i}$ set to zero, so that

$\mathcal{H}_{P}=\sum_{(i,j)\in C_{L}}J_{ij}s_{i}s_{j},$ (9)

with $C_{L}$ denoting available connectivity graph of an D-Wave QPU. The input
classes are:

  * •

RAN1 - Bimodal spin glasses. These inputs are generated by randomly assigning
each coupling a value in $\\{\pm 1\\}$ with equal probability.

  * •

AC3 - Anti-cluster inputs. These inputs are RAN1 problems with couplings
between unit cells multiplied by 3.

  * •

DCL - Deceptive Cluster Loops. The DCL problems Mandrà and Katzgraber (2018)
are a variant of the frustrated cluster loops (FCL) King et al. (2017). Here
we set the parameters of the DCL problems to $(\alpha_{\rm DCL},R_{\rm
DCL},\lambda)=(0.25,1,7)$.

For each size $L$ we generated 50 instances of each RAN1, AC3 and DCL type.

Figure 1: Diversity VS time for three problems on $C_{12}$, or 1152 variables.
The bands show the region between $5^{\rm th}$ and $95^{\rm th}$ percentiles
for the lower bound on diversity. Horizontal black lines show the diversity
target calculated from 50k samples from the QPU. We can readily see the
different behavior between QA and PT+ICM. PT+ICM spends a considerable amount
of time for burn in before settling in good quality basins and starting to
explore good quality and diverse solutions. After that, the diversity of the
samples quickly raises until a plateau whose height depends, as we explain in
the main text, on the number of independent instances of PT+ICM. Adding PT+ICM
instances moves the blue bands up and on the right as more time is spent
burning in (because we consider sequential copies of PT+ICM). Conversely, QA
does not need to burn in and the diversity of the samples returned steadily
increases. From the three problem instances we see that QA is competitive for
all the problems and outperform PT+ICM for AC3 problems. When short times and
lower target diversities are needed, QA clearly outperforms the classical
solver. For longer times, or larger target diversities, and for RAN1 and DCL
inputs, PT+ICM is able to outperform QA.

Figure 2: Comparison of the time-to-diversity (TTD) for parallel tempering
(PT), PT with iso-energetic cluster moves (PT+ICM) and quantum annealing (QA)
for problem built on a Chimera graph with size $L=8$ (top plots) $L=12$
(center plots) and $L=16$ (bottom plots). The instances are ordered according
to increasing TTD. Instances in the grey areas mean that the algorithm reached
the wall-time and it was not possible to calculate the exact TTD. On RAN1
problems we see that although the classical solvers outperforms QA at small
scale, at large scales QA wins over the classical solvers. An even larger
advantage of QA over classical solvers happens for AC3 inputs. For DCL
problems, the situation is somewhat different and PT+ICM has a slight
advantage over QA.

Figure 3: Instance by instance comparison of the Time-To-Diversity (TTD)
between the the D Wave 2000Q solver and PT (upper plots) and PT+ICM (lower
plots), on the three problem instances RAN1, AC3 and DCL on a Chimera size
$L=8$ (left plots) $L=12$ (center plots) and $L=16$ (right plots). Points in
the grey regions mean that the target diversity was not reached within wall
time. A point in the upper left triangle means that, on that particular
instance, QA is winning over the respective classical solver.

###  IV.2 Time Measurements

In order to compare TTD between the solvers, we made the following
assumptions. For QA, we ignore the programming time of the D-Wave 2000Q chip
and consider only the annealing time. For the classical algorithms, we ignore
the time spent for setting the temperatures, and assume single threaded CPU
that can perform a single spin update in 1 ns, so that the elapsed time is
given by measuring the total number of spin updates performed. When referring
to wall-times for classical solvers however, we report actual CPU times on our
machines, taken only by the solver. It is worth noticing that computing
diversity takes much longer than all solver runtimes.

###  IV.3 Experimental Setup

Our experiments consist of four parts: minimum and maximum energy search,
target diversity determination, solvers optimization and performance
comparison. Below we describe each part.

#### Minimum and maximum energy search.

We approximate the ground state $E_{0}$ with the minimum energy $E_{\rm min}$
among a set of up to $2\times 10^{6}$ candidate solutions from the PT+ICM and
QA solvers. For the maximum energy, we used the fact that for bipartite graphs
with zero field $h_{i}=0$, $E_{\rm max}=-E_{0}$.

#### Target Diversity Determination.

As already mentioned, to estimate the target diversity we use the method
suggested in Sect. III.2 and calculate the diversity of 50k samples from the
QPU. The diversity is taken as the maximum of 100 independent calculations
over the 50k samples.

#### Solvers optimization.

We optimized the solvers by a grid search over their hyper-parameter space.
For all solvers we chose a subset of the generated problems, and computed TTD
for each parameter choice. PT and PT+ICM we tuned the exchange rate to set the
temperatures and the number of independent runs. We collect the samples from
every replica in each PT initialization after a fixed number of sweeps. We
then calculate diversity using all samples and compare with the target
diversity. We continue this process until the target diversity is reached. We
found that for both PT and PT+ICM the optimal exchange rate is around $\sim
0.4-0.5$, similarly to the usually adopted exchange rate used for
optimization, while the number of independent runs is $\sim 25$ for PT and
$\sim 10$ for PT+ICM. For the latter however, we note that the ICM moves
include an extra factor of 2 in the number of replicas.

For QA we optimized the annealing time as well as the number of spin-reversal
transformations. We found that the lowest possible annealing time of 1 $\mu$s
minimized the ${\rm TTD}$, while increasing the number of spin-reversal
transformations helps to increase the diversity.

Details of the optimization are shown in App. B.

###  IV.4 Solvers Behavior

Before presenting the results of our analysis, we briefly discuss behaviors of
the QPU and PT+ICM on the three types of problem (PT behaves similarly as
PT+ICM and as we show in the next section, performs slightly worse than
PT+ICM). In Fig. 1 we plot the diversity reached by a solver over time for
three problems - one for each input class - on the $C_{12}$ graph. The bands
show the $5^{\rm th}$ and $95^{\rm th}$ percentiles of the distribution of
diversity, while the black solid lines are the diversity threshold set by the
50k QPU samples. Here we can see the fundamental difference between the two
types of solvers. PT+ICM has low diversity during the burn-in time, then
quickly increases when each instance of PT access low energy configurations
and eventually saturates, possibly because it remains stuck in a limited
region. On the contrary, QA is able to produce a set of diverse solutions very
quickly and the diversity of the sample set increases steadily with time. More
importantly, during the times explored in this analysis, QA does not saturate
meaning that the QPU does not get stuck in a certain region in the solution
space.

From Fig. 1 it is also possible to qualitatively inspect what happens when the
target diversity is changed. By increasing the number of samples of the D Wave
solver, we effectively raise the horizontal lines for each problem, the exact
amount depending on the steepness of the orange bands. In all three problems,
increasing the target diversity above the saturation level, the classical
solver with the current (optimized) setup will not be able to reach the target
diversity obtained by the QPU within a reasonable time. This however may be
fixed by another choice of hyper-parameters, for example by increasing the
number of copies per temperature, but that will also increase TTD as the
classical solver will need more total time to converge.333With this we mean
that although each copy spends approximately the same time to converge, having
more copies considerably increase the cumulative time spent. Moreover, one
also needs to take into account the time required to find a new set of hyper-
parameters for each setting.

For the specific RAN1 instance showed here, an increase of target diversity of
about 15% - obtained using approximately 100k samples from the QPU - would
degrade the relative advantage of the QPU over PT+ICM. For the AC3 problem
shown in the center panel, QPU has an advantage at all the times shown.
Increasing or decreasing the samples used to calculate the target diversity
would increase the advantage. For the DCL problem, increasing the diversity
threshold would make PT+ICM win over QA until the plateau at $t\gtrsim 3\times
10^{5}$ $\mu$s. Beyond that limit, as already mentioned, one should re-
optimize the hyper-parameters of the solvers and investigate the solution
again.

###  IV.5 Performance Comparison

Here we show the performance comparison of the solvers. For each solver and
for each problem instance we performed 100 independent experiments, each of
them consisting of measuring diversity over computation time, as described in
Sect. IV.2. For the classical solvers we imposed a CPU wall time of 5 minutes,
while for the QPU we imposed a maximum number of analyzed samples of $10^{6}$.
The main results are shown in Figs. 2 and 3. In the former we show the results
for the three solvers, ordered by increasing TTD from left to right. The left
panels show the results for the RAN1 input instances, the center are the AC3
instances, and the right panels are the DCL inputs. Moving vertically we
change the size of the Chimera graph. For the AC3 instances we see that the
QPU performs better than the other classical solvers, with QPU relative
performances improving with larger sizes. This can also be seen in Fig. 3,
where we show instance by instance comparison between the solvers. On the top
panels we show the comparison between PT and QA, while on the lower panels we
compare the QPU with PT+ICM. Moving horizontally we increase the size of the
Chimera graph.

In both figures, we see that for AC3 instances (yellow triangles) the QPU is
faster than the two classical solvers and how for increasing sizes the
instance by instance comparisons “move up” in the plots. Likewise for RAN1
instances, the QPU has better performances, although on $C_{8}$ the classical
solvers have comparable, if not better performances. On C8 in fact, as can be
seen in the lower left panel on Fig. 2, PT+ICM is faster than QA. We note that
for RAN1 and AC3, the TTD for QPU is almost constant among the instances. This
is probably due to the fact that the target diversity is set by the QPU. For
DCL problems the behaviour is different and the performance of the QPU is
comparable with that of PT and PT+ICM and does not get better with increasing
size. This could be related to the fact that DCL parameters require more
precision.

###  IV.6 Decreasing the approximation ratio

Figure 4: Comparison of TTD with approximation ratio $\alpha=0.005$ for PT+ICM
and QA for problems on $C_{8}$ (upper panels), $C_{1}2$ (mid panels) and
$C_{1}6$ (lower panels). From left to right the inputs are RAN1, AC3 and DCL.
Here the situation is similar to Fig. 2: for RAN1 inputs, QA performance
relative to PT+ICM improves as the scale of the problem is increased. For AC3
QA outperforms the classical solver at all scales and for DCL problems,
interestingly, the performances are on par at all scales.

Here we briefly discuss what happens when a lower approximation ratio is
required. For lower $\alpha$, TTD will be dominated by the TTT metric King et
al. (2015b). We focused here on $\alpha=0.005$ and to limit computational
time, we selected only a small batch of the problems above. Also, because
PT+ICM performed consistently better than PT, we only considered the former.
We re-optimzed the solvers and found that for PT+ICM the optimal exchange rate
dropped to $0.2$ while the number of copies per temperature remained the same,
10. However, in this case the optimal region in the parameter space is larger,
as can be seen in Fig. 8. For the QPU the optimal parameters did not change.

The results for this approximation ratio are shown in Fig. 4 and they are
similar to the case of $\alpha=0.01$. For AC3 problems, TTD for QA is
consistently smaller than for PT+ICM. In the case of RAN1 problems, PT+ICM
performs better for smaller problems, but TTD scales in favour of the QPU for
larger sizes. For DCL inputs, the performances between the two solvers are
comparable for all sizes. Due to calibration errors, noise and finite
temperature, we expect the performance of the QPU do degrade as we decrease
the approximation ratio to lower values. Although $\alpha=0.005$ is already a
good target for many applications, future generations of the QPU are likely to
have less calibration errors and lower noise, making smaller approximation
ratios achievable.

##  V Conclusions

In this work, we have used a new metric to assess the diversity of the samples
returned by quantum annealers. In principle, many definitions of diversity are
possible, from clustering-inspired metrics to other metrics like average
Hamming distance between a sequence of solutions. Here, the diversity is the
maximum independent set of the graph of samples below an approximation ratio
$\alpha$, where the edges are defined only between samples within a Hamming
radius $R$. We believe that such a diversity metric can have a practical
relevance in the context of hybrid genetic algorithms and quality diversity
optimization Pugh et al. (2016). Although computing the metric is #P-hard, we
identified an algorithm that can quickly compute its lower bound, i.e. it
scales linearly with the number of solutions. Similarly to the Time-To-
Solution (TTS) metric, we have adopted the notion of Time-To-Diversity (TTD)
as a metric to assess the speed of a heuristic to obtain a target diversity of
a problem. We then assessed the performance of the solvers to reach the
reference diversity for the case $\alpha=0.01$ and $R=0.25$. We observed that
for RAN1 and AC3 problems, the D Wave quantum annealer shows a performance
advantage over the classical solvers, with the advantage increasing with
problem size. For the DCL inputs, the QPU performances are comparable with
PT+ICM at all sizes. We also considered a case in which the approximation
ratio is reduced to $\alpha=0.005$ and found that the results were similar to
the case of $\alpha=0.01$. This suggests that, for applications were large
diversity of samples is required, a portfolio solver that combines quantum and
classical heuristics will win over other solvers.

###### Acknowledgements.

This work was supported by Mitacs through the Mitacs Accelerate program. The
authors would like to thank Jeremy Hilton, Emile Hoskinson, Andrew King, Eric
Ladizinsky, Jack Raymond and Trevor Lanting for useful discussions.

## References

  * Mohseni et al. (2021) M. Mohseni, M. M. Rams, S. V. Isakov, D. Eppens, S. Pielawa, J. Strumpfer, S. Boixo, and H. Neven, arXiv e-prints arXiv:2110.10560 (2021), eprint 2110.10560.
  * Johnson et al. (2011) M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, et al., Nature 473 (2011), ISSN 0028-0836.
  * Zhu et al. (2015) Z. Zhu, A. J. Ochoa, and H. G. Katzgraber, Phys. Rev. Lett. 115, 077201 (2015), eprint 1501.05630.
  * Zhu et al. (2016) Z. Zhu, C. Fang, and H. G. Katzgraber, arXiv e-prints arXiv:1605.09399 (2016), eprint 1605.09399.
  * Hen et al. (2015) I. Hen, J. Job, T. Albash, T. F. Rønnow, M. Troyer, and D. A. Lidar, Physical Review A 92, 042325 (2015), eprint 1502.01663.
  * King et al. (2015a) A. D. King, T. Lanting, and R. Harris, arXiv e-prints arXiv:1502.02098 (2015a), eprint 1502.02098.
  * King et al. (2017) J. King, S. Yarkoni, J. Raymond, I. Ozfidan, A. D. King, M. M. Nevisi, J. P. Hilton, and C. C. McGeoch, arXiv e-prints arXiv:1701.04579 (2017), eprint 1701.04579.
  * Albash and Lidar (2018) T. Albash and D. A. Lidar, Physical Review X 8, 031016 (2018), eprint 1705.07452.
  * Mandrà and Katzgraber (2018) S. Mandrà and H. G. Katzgraber, Quantum Science and Technology 3, 04LT01 (2018), eprint 1711.01368.
  * Dickson et al. (2013) N. G. Dickson, M. W. Johnson, M. H. Amin, R. Harris, F. Altomare, A. J. Berkley, P. Bunyk, J. Cai, E. M. Chapple, P. Chavez, et al., Nature communications 4 (2013), ISSN 2041-1723.
  * Boixo et al. (2014) S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Nature Physics 10 (2014), ISSN 1745-2473.
  * Katzgraber et al. (2015) H. G. Katzgraber, F. Hamze, Z. Zhu, A. J. Ochoa, and H. Munoz-Bauza, Physical Review X 5 (2015), ISSN 2160-3308.
  * Heim et al. (2015) B. Heim, T. F. Rønnow, S. V. Isakov, and M. Troyer, Science (New York, N.Y.) 348 (2015), ISSN 1095-9203.
  * Katzgraber et al. (2014) H. G. Katzgraber, F. Hamze, and R. S. Andrist, Physical Review X 4, 021008 (2014), eprint 1401.1546.
  * Katzgraber et al. (2015) H. G. Katzgraber, F. Hamze, Z. Zhu, A. J. Ochoa, and H. Munoz-Bauza, Physical Review X 5, 031026 (2015), eprint 1505.01545.
  * Heim et al. (2015) B. Heim, T. F. Rønnow, S. V. Isakov, and M. Troyer, Science 348, 215 (2015), eprint 1411.5693.
  * Albash et al. (2015) T. Albash, T. F. Rønnow, M. Troyer, and D. A. Lidar, European Physical Journal Special Topics 224 (2015).
  * Martin-Mayor and Hen (2015) V. Martin-Mayor and I. Hen, Scientific Reports 5, 15324 (2015), eprint 1502.02494.
  * Marshall et al. (2016) J. Marshall, V. Martin-Mayor, and I. Hen, Physical Review A 94, 012320 (2016), eprint 1605.03607.
  * Denchev et al. (2016) V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy, J. Martinis, and H. Neven, Physical Review X 6, 031015 (2016), eprint 1512.02206.
  * King et al. (2021) A. D. King, J. Raymond, T. Lanting, S. V. Isakov, M. Mohseni, G. Poulin-Lamarre, S. Ejtemaee, W. Bernoudy, I. Ozfidan, A. Y. Smirnov, et al., Nature Communications 12, 1113 (2021).
  * King et al. (2015b) J. King, S. Yarkoni, M. M. Nevisi, J. P. Hilton, and C. C. McGeoch, arXiv e-prints arXiv:1508.05087 (2015b), eprint 1508.05087.
  * Matsuda et al. (2009) Y. Matsuda, H. Nishimori, and H. G. Katzgraber, in _Journal of Physics Conference Series_ (2009), vol. 143 of _Journal of Physics Conference Series_ , p. 012003, eprint 0808.0365.
  * Mandrà et al. (2017) S. Mandrà, Z. Zhu, and H. G. Katzgraber, Phys. Rev. Lett. 118, 070502 (2017), eprint 1606.07146.
  * Könz et al. (2018) M. S. Könz, G. Mazzola, A. J. Ochoa, H. G. Katzgraber, and M. Troyer, arXiv e-prints arXiv:1806.06081 (2018), eprint 1806.06081.
  * Zhang et al. (2017) B. H. Zhang, G. Wagenbreth, V. Martin-Mayor, and I. Hen, Scientific Reports 7, 1044 (2017), eprint 1701.01524.
  * Zhu et al. (2019) Z. Zhu, A. J. Ochoa, and H. G. Katzgraber, Phys. Rev. E 99, 063314 (2019), eprint 1903.07600.
  * King et al. (2019) J. King, S. Yarkoni, J. Raymond, I. Ozfidan, A. D. King, M. M. Nevisi, J. P. Hilton, and C. C. McGeoch, Journal of the Physical Society of Japan 88, 061007 (2019), eprint https://doi.org/10.7566/JPSJ.88.061007, URL https://doi.org/10.7566/JPSJ.88.061007.
  * Pugh et al. (2016) J. K. Pugh, L. B. Soros, and K. O. Stanley, Frontiers in Robotics and AI 3, 40 (2016), ISSN 2296-9144, URL https://www.frontiersin.org/article/10.3389/frobt.2016.00040.
  * Ulrich et al. (2010) T. Ulrich, J. Bader, and E. Zitzler, in _Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation_ (Association for Computing Machinery, New York, NY, USA, 2010), GECCO ’10, p. 455–462, ISBN 9781450300728, URL https://doi.org/10.1145/1830483.1830569.
  * Ulrich and Thiele (2011) T. Ulrich and L. Thiele, in _Proceedings of the 13th Annual Conference on Genetic and Evolutionary Computation_ (Association for Computing Machinery, New York, NY, USA, 2011), GECCO ’11, p. 641–648, ISBN 9781450305570, URL https://doi.org/10.1145/2001576.2001665.
  * Gao et al. (2015) W. Gao, S. Nallaperuma, and F. Neumann, arXiv e-prints arXiv:1510.08568 (2015), eprint 1510.08568.
  * Alexander et al. (2017) B. Alexander, J. Kortman, and A. Neumann, in _Proceedings of the Genetic and Evolutionary Computation Conference_ (Association for Computing Machinery, New York, NY, USA, 2017), GECCO ’17, p. 171–178, ISBN 9781450349208, URL https://doi.org/10.1145/3071178.3071342.
  * Neumann et al. (2018) A. Neumann, W. Gao, C. Doerr, F. Neumann, and M. Wagner, arXiv e-prints arXiv:1802.05448 (2018), eprint 1802.05448.
  * Mohseni et al. (2021) M. Mohseni et al., in prep. (2021).
  * Bunyk et al. (2014) P. I. Bunyk, E. M. Hoskinson, M. W. Johnson, E. Tolkacheva, F. Altomare, A. J. Berkley, R. Harris, J. P. Hilton, T. Lanting, A. J. Przybysz, et al., IEEE Transactions on Applied Superconductivity 24, 2318294 (2014), eprint 1401.5504.
  * Lucas (2014) A. Lucas, Frontiers in Physics 2, 5 (2014), eprint 1302.5843.
  * Boothby et al. (2016) T. Boothby, A. D. King, and A. Roy, Quantum Information Processing 15, 495 (2016), eprint 1507.04774.
  * Cai et al. (2014) J. Cai, W. G. Macready, and A. Roy, arXiv e-prints arXiv:1406.2741 (2014), eprint 1406.2741.
  * Ehrenfest (1916) P. Ehrenfest, Annalen der Physik 356, 327 (1916), URL https://onlinelibrary.wiley.com/doi/abs/10.1002/andp.19163561905.
  * Born and Fock (1928) M. Born and V. Fock, Zeitschrift für Physik 51, 165 (1928).
  * Saket (2013) R. Saket, arXiv e-prints arXiv:1306.6943 (2013), eprint 1306.6943.
  * Kosowski and Manuszewski (2004) A. Kosowski and K. Manuszewski (2004).
  * Matula and Beck (1983) D. W. Matula and L. L. Beck, J. ACM 30, 417 (1983), ISSN 0004-5411, URL http://doi.acm.org/10.1145/2402.322385.
  * Hukushima (1999) K. Hukushima, Phys. Rev. E 60, 3606 (1999), eprint cond-mat/9903391.

##  Appendix A Algorithms for approximated diversity

Here we provide the details of the algorithms used to compute lower and upper
bounds on diversity.

###  A.1 Lower bound: left-neighbors-aggregation algorithm

We begin with the agglomerative algorithm used to compute the lower bound on
the diversity. Given the samples returned by a solver, we only keep the unique
fit samples in order of appearance. This defines a _labelled_ graph
$G=(V,L,E)$ where $L:V\to\\{1,\dots,m\\}$ is the _labelling_ that associates
each node $v\in V$ with an integer number, $L(v)\in\\{1,\dots,m\\}$ \- in our
case, given by the order in which the sampler is returned by the solver. Given
a vertex $v$, we also define its neighbours set N$(v)$ as the set of nodes
connected to $v$ through an edge $e\in E$. We also define the left (right)
neighbours LN$(v)$ (RN$(u)$) as the subset of the neighbours $u$ of $v$ such
that $L(u)<L(v)$ ($L(u)>L(v)$). The algorithm is illustrated in Alg. 1.

It is clear that the set $\mathcal{S}_{\alpha}^{*}$ constructed by the LNA
algorithm depends on the specific labelling function, and that it can also be
a maximum independent set. The latter property is easy to show. Assuming that
the graph has a maximum independent set of $\alpha$ elements, one can just
pick a labelling $L$ for which the first $\alpha$ elements are the elements of
the maximum independent set, and then the other nodes of the graph. The
advantage of using this algorithm to find a lower bound on diversity is that
its time-complexity is $\mathcal{O}(m\cdot D^{\rm
low}_{R\mathcal{S}_{\alpha}})$ at maximum, making it suitable for our
purposes.

To improve the lower bound on diversity given by this algorithm, we could re-
shuffle the collected samples.

Algorithm 1 Left-Neighbors-Aggregation

INPUT: $\mathcal{S}_{\alpha}$: set of fit, unique samples.
$|\mathcal{S}_{\alpha}|=m$. Labelling function
$L:\mathcal{S}_{\alpha}\to{1,\dots,m}$

Set $\mathcal{S}^{*}_{\alpha}=\left\\{L^{-1}(1)\right\\}$

for all $i=2,\dots,n$ do

Set $v=L^{-1}(i)$

if ${\rm LN}(v)\cap\mathcal{S}_{\alpha}^{*}=\varnothing$ then

$\mathcal{S}_{\alpha}^{*}=\mathcal{S}_{\alpha}^{*}\cap{v}$

end if

end for

OUTPUT: $\mathcal{S}_{\alpha}^{*}$.

###  A.2 Upper bound: greedy coloring algorithm

Here we describe the _greedy-coloring_ algorithm Kosowski and Manuszewski
(2004) used to provide the upper bound on diversity. As we did with the
agglomerative algorithm in App. A.1, we keep all the samples returned by a
solver that have energy within the threshold $E_{\rm th}$ in order of
appearance. This defines the complement graph $\bar{G}=(V,L,\bar{E})$ where
$L$ is the labelling function defined in App. A.1, and $\bar{E}$ is the set of
edges defined between vertices (samples) that have Hamming distance greater
than $nR$. The greedy-coloring algorithm then follows the ordering given by
$L$ and is summarized in Alg. 2.

Algorithm 2 Greedy - coloring

INPUT: $\mathcal{S}_{\alpha}$: set of fit, unique samples.
$|\mathcal{S}_{\alpha}|=m$. Labelling function
$L:\mathcal{S}_{\alpha}\to{1,\dots,m}$

Set $C_{1}={L^{-1}(1)}$, $c=1$

for all $i=2,\dots,n$ do

Set $v=L^{-1}(i)$

Set $j=1$

vertex $v$ is not assigned

while (vertex $v$ is not assigned) $j\leq c$ do

if $C_{j}\cap N(v)=\varnothing$ then

$C_{j}=C_{j}\cup\\{v\\}$

$v$ is assigned

else

j = j+1

end if

end while

if $v$ is not assigned then

$C_{c+1}=\\{v\\}$

$c=c+1$

end if

end for

OUTPUT: number of colors $c$.

Like the agglomerative algorithm in App. A.1, the greedy-coloring depends on
the labelling $L$. Different strategies to choose the labelling can lead to
considerable improvement of performances.

###  A.3 Algorithms performances

Here we compare the performance of the LNA and greedy coloring algorithms
against another tight lower bound on diversity, i.e. by mapping the maximum
independent set problem to an Ising Hamiltonian Lucas (2014) and solving with
via simulated annealing (SA). In Fig. 5 we show the comparison of the
algorithms for one problem in each class analyzed in this paper. The mapping
to Ising Hamiltonian is denoted by “exact” meaning that we are solving the
exact maximum independent set problem with the SA heuristic. To improve the
accuracy of the LNA algorithm, we reshuffle the samples 100 times. We can see
how with this setting - used also in our analysis in Sect. IV the lower bound
is remarkably close to the values obtained by solving the maximum independent
set problem. Conversely, the upper bound is not a tight bound and diverges
from the exact value rather quickly. Importantly, the computation time for the
LNA method has a better scaling compared to the other two methods. This
motivates us to adopt the lower bound quantity for our analysis.

Figure 5: Comparison between the approximated diversity algorithms LNA and
greedy coloring and a tight bound on diversity. The upper panels show the
diversity computed by algorithms as a function of samples analyzed for three
instances in each class RAN1, AC3 and DCL. The lower panels show the CPU time
spent on the computations.

##  Appendix B Solvers optimization

In this Appendix we describe the solver optimization.

###  B.1 Parallel Tempering

Figure 6: Configuration of PT (left side) and PT+ICM (right side) used in this
work. Black vertical arrows denote replica exchange due to thermal
fluctuations, while red horizontal arrows denote the Houdayer iso-energetic
cluster moves, performed between replicas at the same temperature.

Figure 7: Medians of TTD ($\alpha=0.01$, $R=0.25$) for a selection of problems
on C12, for the hyper-parameters of PT (left panel) and PT+ICM (right panel).
The “x” mark shows the optimal configuration of the hyper-parameters. Color
coded in logarithmic scale with blue standing for largest TTD and yellow for
smallest TTD.  Figure 8: Medians of TTD ($\alpha=0.005$, $R=0.25$) for a
selection of problems on C12, for the hyper-parameters of PT+ICM. The “x” mark
shows the optimal configuration of the hyper-parameters. Color coded in
logarithmic scale with blue standing for largest TTD and yellow for smallest
TTD. For this approximation ratio we did not consider the PT solver. With
respect to the case in Fig. 7 the degeneracy between the two parameters is
more pronounced.

For PT and PT+ICM we chose to optimize the commonly used target rate (the
exchange rate between neighboring replicas) and the number of copies per
temperature. For the former we adjusted the number and the spacing of the
temperatures using the method of Hukushima (1999). In our experiments we found
beneficial to collect samples from multiple copies at the same temperature,
see Fig. 6. These can be seen as independent PT solvers running in parallel.
However, there is a trade-off in the “number of copies” of PT running as while
increasing them helps increasing the diversity of samples, the overhead of
having multiple PTs considerably increases TTD. The results of PT
optimizations are shown in Fig. 7. We can see a degeneracy between target rate
- proportional to the number of replicas - and number of copies of PT. This is
expected as more copies of PT consist in an effective increase in number of
replicas. Interestingly however, having only 1 copy of PT with larger target
rate result in consistent larger TTD.

###  B.2 Quantum Annealing

Figure 9: Optimization of quantum annealing. The colors show the total
diversity after 50k samples: blue standing for low diversity and yellow for
large diversity. The diversity is computed for an approximation ratio
$\alpha=0.01$ and a radius factor $R=0.25$.

For quantum annealing we optimized the annealing time and the number of spin
reversal transformations. We chose a subset of the inputs and optimized
overall diversity with approximation ratio $\alpha=0.01$ and radius factor
$R=0.25$. Spin-reversal transformations are local transformations
$\alpha=[\alpha_{1},\dots,\alpha_{n}]\in\\{\pm 1\\}^{\times n}$ that flip the
spin $i$ if $\alpha_{i}=-1$, accompanied by the Hamiltonian transformations

$h_{i}\to\alpha_{i}h_{i},\quad J_{ij}\to\alpha_{i}\alpha_{j}J_{ij}.$ (10)

Although the Ising Hamiltonian (1) does not change under these
transformations, they can reduce systematic and analog errors on the QPU. We
compute the total diversity of 50k samples obtained by the QPU with hyper-
parameters on the grid $t_{a}\in[1\,\mu{\rm s},5\,\mu{\rm s},10\,\mu{\rm s}]$
and number of spin reversal transformations (for each 10k samples batch) in
[0, 1, 10, 100]. The results are shown in Fig. 9.


#  The eigenvectors-eigenvalues identity and Sun’s conjectures on determinants
and permanents

Xuejun Guo  Xuejun Guo
Department of Mathematics
Nanjing University
Nanjing 210093, China<EMAIL_ADDRESS>,  Xin Li  Xin Li
Department of Mathematics
Nanjing University
Nanjing 210093, China<EMAIL_ADDRESS>,  Zhengyu Tao  Zhengyu Tao
Department of Mathematics
Nanjing University
Nanjing 210093, China<EMAIL_ADDRESS>and  Tao Wei  Tao Wei
Department of Mathematics
Nanjing University
Nanjing 210093, China<EMAIL_ADDRESS>
###### Abstract.

In this paper, we prove a conjecture raised by Zhi-Wei Sun in 2018 by the
eigenvectors-eigenvalues identity found by Denton, Parke, Tao and X. Zhang in
2019.

The authors are supported by National Nature Science Foundation of China (Nos.
11971226, 11631009).

Keywords: eigenvectors-eigenvalues identity, determinant

2020 Mathematics Subject Classification:Primary 11C20, 15A15; Secondary 05A19,
11A07, 33B10

##  1\. Introduction

In 2019, Peter B. Denton, Stephen J. Parke, Terence Tao and Xining Zhang [8]
found the legendary eigenvector-eigenvalue identity and showed how this
identity can also be used to extract the relative phases between the
components of any given eigenvector. In this paper, we will prove several
conjectural identities of Zhi-Wei Sun on determinants and permanents of
matrices involving roots of unity. We will extract the eigenvalues of the
minors $M_{j}(1\leq j\leq n)$ of the Hermitian matrix $A$ from the
eigenvectors of $A$ by virtue of the eigenvector-eigenvalue identity. So it
can been seen as extracting eigenvalues from eigenvectors.

Zhi-Wei Sun began a systematical research on permanents and determinants of
matrices in number theory in 2018. He found many identities with rich
arithmetic meanings and also raised many open problems. His results are in
[3], [4], [5], [6]. One of Sun’s conjectures is the following one which was
firstly raised in 2018.

######  Conjecture 1.1 (Zhi-Wei Sun, Conjecture 4.3 of [6]).

Let $n>1$ be a integer and $\zeta$ a primitive $n$-th root of unity.

  1. (1)

If $n$ is a even number, then

$\sum_{\tau\in
D(n)}\prod_{j=1}^{n}\dfrac{1}{1-\zeta^{j-\tau(j)}}=\dfrac{\left(\left(n-1\right)!!\right)^{2}}{2^{n}}=\dfrac{n!}{4^{n}}\binom{n}{n/2},$
(1.1)

  2. (2)

if $n$ is an odd number, then

$\sum_{\tau\in
D(n-1)}\prod_{j=1}^{n-1}\dfrac{1}{1-\zeta^{j-\tau(j)}}=\dfrac{1}{n}\left(\dfrac{n-1}{2}!\right)^{2},$
(1.2)

and

$\sum_{\tau\in
D(n-1)}\operatorname{sign}(\tau)\prod_{j=1}^{n-1}\dfrac{1}{1-\zeta^{j-\tau(j)}}=\frac{(-1)^{\frac{n-1}{2}}}{n}\left(\frac{n-1}{2}!\right)^{2},$
(1.3)

where $D(n-1)$ is the set of all derangements $\tau$ of indices
$j=1,\ldots,n-1$ such that $\tau(j)\neq j$ for all $j=1,\ldots,n-1$.

##  2\. Eigenvalues from eigenvectors

In this section, we will prove the equation (1.3). The left hand side of (1.3)
is the determinant of an $(n-1)\times(n-1)$ sub-matrix of the $n\times n$
matrix $A=(a_{jk})$ such that
$a_{jk}=\left(1-\delta_{jk}\right)\left\\{1+i\cot\frac{(j-k)\pi}{n}\right\\}$,
where $\delta_{jk}=1$ if $j=k$; and $0$ if $j\neq k$. The following theorem
gives the eigenvalues and eigenvectors of $A$.

######  Theorem 2.1 (F. Calogero and A.M. Perelomov, Theorem 1 of [1]).

The off-diagonal hermitian matrix of order $n$ whose elements are defined by
the formula

$a_{jk}=\left(1-\delta_{jk}\right)\left\\{1+i\cot\frac{(j-k)\pi}{n}\right\\}$

has the integer eigenvalues

$\lambda_{i}=2i-n-1,\quad i=1,2,\ldots,n$

and the corresponding eigenvectors $v_{i}$ with norm $1$ have components

$v_{i,j}=\exp\left(-\frac{2\pi{\rm{i}}ij}{n}\right)/\sqrt{n},\quad
j=1,2,\ldots,n.$

If $A$ is an $n\times n$ Hermitian matrix, then we denote its $n$ real
eigenvalues by $\lambda_{1}(A),\ldots,\lambda_{n}(A)$. We can find an
orthonormal basis $v_{1},v_{2},\cdots,v_{n}\in\mathbb{C}^{n}$ such that
$v_{i}$ is associated to $\lambda_{i}(A),1\leq i\leq n$. For any
$i,j=1,\ldots,n$, let $v_{i,j}$ denote the $j^{\text{th }}$ component of
$v_{i}.$ If $1\leq j\leq n$, let $M_{j}$ denote the $(n-1)\times(n-1)$ minor
formed from $A$ by deleting the $j^{\text{th }}$ row and column from $A$. This
is also a Hermitian matrix, and thus has $n-1$ real eigenvalues
$\lambda_{1}\left(M_{j}\right),\ldots,\lambda_{n-1}\left(M_{j}\right)$.

######  Theorem 2.2 (Eigenvector-eigenvalue identity, Theorem 1 of [8]).

With the notation as above, we have

$\left|v_{i,j}\right|^{2}\prod_{k=1;k\neq
i}^{n}\left(\lambda_{i}(A)-\lambda_{k}(A)\right)=\prod_{k=1}^{n-1}\left(\lambda_{i}(A)-\lambda_{k}\left(M_{j}\right)\right).$
(2.1)

Proof of Equation (1.3):

Note that

$\displaystyle\frac{1}{1-\zeta^{j-k}}$
$\displaystyle=\frac{1}{2}\left(1+\frac{1+\zeta^{j-k}}{1-\zeta^{j-k}}\right)$
$\displaystyle=\frac{1}{2}\left(1+{\rm{i}}\cot\frac{\pi(j-k)}{n}\right).$

Let $A=(A_{jk})$ be the $n\times n$ Hermitian matrix defined in Theorem 2.1.
By Theorem 2.1, the eigenvalues of $A$ are

$\lambda_{1}=n-1,\lambda_{2}=n-3,\cdots,\lambda_{i}=n+1-2i,\cdots,\lambda_{n}=1-n.$
(2.2)

Let $i=\frac{n+1}{2}$ and $j=n$. Then $\lambda_{i}=0$. By Theorem 2.2,

$\left|v_{i,n}\right|^{2}\prod_{k=1;k\neq
i}^{n}\left(-\lambda_{k}(A)\right)=\prod_{k=1}^{n-1}\left(-\lambda_{k}\left(M_{n}\right)\right).$

By Theorem 2.1, $\left|v_{s,n}\right|^{2}=\frac{1}{n}$. Hence

$|M_{n}|={\prod_{k=1}^{n-1}\lambda_{k}\left(M_{n}\right)}=\frac{(-1)^{\frac{n-1}{2}}((n-1)!!)^{2}}{n}.$

Hence

$\displaystyle\sum_{\tau\in
D(n-1)}\operatorname{sign}(\tau)\prod_{j=1}^{n-1}\frac{1}{1-\zeta^{j-\tau(j)}}$
$\displaystyle=2^{1-n}|M_{n}|=\frac{2^{1-n}(-1)^{\frac{n-1}{2}}((n-1)!!)^{2}}{n}$
$\displaystyle=\frac{(-1)^{\frac{n-1}{2}}}{n}\left(\frac{n-1}{2}!\right)^{2}.$

Hence Equation (1.3) is proved. ∎

Han Wang and Zhi-Wei Sun prove [7] another conjecture involving derangements
and roots of unity by similar arguements. My colleague Keqin Liu suggested a
different strategy to prove Conjecture 1.1. Although the eigenvalues of $A$
are rational integers as in Theorem 2.1, the Numerical computation shows that
there is no simple pattern of the eigenvalues of $M_{n}$ in the proof of
equation (1.3) . However the eigenvalues of $M_{n}B$ are very simple for
carefully chosen diagonal matrices. Let $B_{s}$ the diagonal matrix whose
diagonal entries are given as $1-\zeta^{is},\quad 1\leq i\leq n-1,$ where
$s\in\left\\{-\frac{n-1}{2},-\frac{n-3}{2},\cdots,-1,1,2,\cdots,\frac{n-1}{2}\right\\}$.
Let $M_{n,s}=M_{n}B_{s}$. Then (1.3) is equivalent to

$\operatorname{det}\left(M_{n,1}\right)=(-1)^{\frac{n-1}{2}}\left(\frac{n-1}{2}!\right)^{2}.$
(2.3)

Liu found that the eigenvalues of $M_{n,1}$ for are

$-\frac{n-1}{2},\cdots,-1,1,\cdots,\frac{n-1}{2}$

by numerical computation. He discusses this strategy in [2].

Let $f(x)=\prod_{k=1}^{n-1}\left(x-\lambda_{k}\left(M_{j}\right)\right)$. Then
by (2.1) and (2.2),

$f(\lambda_{i}(A))=f(n+1-2i)=\frac{1}{2n}\prod_{k=1;k\neq
i}^{n}\left(k-i\right),\ 1\leq i\leq n.$ (2.4)

Hence $f(x)$ is completely determined by the Lagrange interpolation.

##  3\. The permanents of matrices

In this section, we will prove equation (1.1) and equation (1.2).

######  Lemma 3.1.

Let $n$ be an integer greater than $2$, and $x_{1},\cdots,x_{n}$ pairwise
distinct complex numbers. Let $L(n)$ be the set of all elements $\tau\in S(n)$
such that $\tau$ is an $n$-cycle. Then

$\sum_{\tau\in L(n)}\prod_{j=1}^{n}\frac{1}{x_{\tau(j)}-x_{j}}=0.$

###### Proof.

For $\tau_{1}=(1\ a_{2}\ a_{3}\ \cdots\ a_{n}),\tau_{2}=(1\ b_{2}\ b_{3}\
\cdots\ b_{n})\in L(n)$, define $\tau_{1}\sim\tau_{2}$ if and only if $(a_{2}\
a_{3}\ \cdots\ a_{n})=(b_{2}\ b_{3}\ \cdots\ b_{n})\in L(n-1)$. One can easily
see that $\sim$ is an equivalence relation on $L(n)$. Obviously, there are
$n-1$ elements in each equivalence class. Take an arbitrary equivalence class

$\displaystyle\sigma_{1}$ $\displaystyle=(1\ a_{2}\ a_{3}\ \cdots\ a_{n-1}\
a_{n}),$ $\displaystyle\sigma_{2}$ $\displaystyle=(1\ a_{3}\ a_{4}\ \cdots\
a_{n}\ a_{2}),$ $\displaystyle\ \ \ \ \ \vdots$ $\displaystyle\sigma_{n-1}$
$\displaystyle=(1\ a_{n}\ a_{2}\ \cdots\ a_{n-2}\ a_{n-1}).$

Then we just need to show
$\displaystyle\sum_{i=1}^{n-1}\prod_{j=1}^{n}\frac{1}{x_{\sigma_{i}(j)}-x_{j}}=0$.
Let
$\displaystyle\lambda_{i}=\prod_{j=1}^{n}\frac{1}{x_{\sigma_{i}(j)}-x_{j}}$,
obviously, $\alpha_{i}\neq 0$ and

$\frac{\lambda_{i}}{\lambda_{1}}=\left\\{\begin{array}[]{ll}\frac{(x_{1}-x_{a_{2}})(x_{1}-x_{a_{n}})}{x_{a_{2}}-x_{a_{n}}}\left(\frac{1}{x_{1}-x_{a_{i+1}}}-\frac{1}{x_{1}-x_{a_{i}}}\right),&i\neq
1,\\\ 1,&i=1.\end{array}\right.$

Thus
$\displaystyle\sum_{i=1}^{n-1}\frac{\lambda_{i}}{\lambda_{1}}=1+\frac{(x_{1}-x_{a_{2}})(x_{1}-x_{a_{n}})}{x_{a_{2}}-x_{a_{n}}}\left(\frac{1}{x_{1}-x_{a_{n}}}-\frac{1}{x_{1}-x_{a_{2}}}\right)=0$.
Hence we have $\displaystyle\sum_{i=1}^{n-1}\alpha_{i}=0.$ This complete the
proof of lemma 3.1. ∎

######  Theorem 3.2.

Let $n$ be a positive integer, and $x$ a complex number such that $x^{k}\neq
1$ for any positive integer $k<n$. Assume

$A=(a_{ij})_{n\times
n}=\displaystyle\left((1-\delta_{ij})\left(1-\frac{1}{1-x^{i-j}}\right)\right)_{n\times
n},$

and

$B=(b_{ij})\in M_{(n-k)\times(n-k)}(\mathbb{C})$

be the sub-matrix obtained from $A$ by deleting the $s_{1}th,\cdots,s_{k}th$
columns and the $s_{1}th,\cdots,s_{k}th$ rows, where $k=0$ or $2\leqslant
k<n,\ 1\leqslant s_{1}<s_{2}<\cdots<s_{k}\leqslant n$.

  1. (1)

If $n-k$ is an odd number, then

$\sum_{\begin{subarray}{c}\tau\in D(n-k)\\\
\mathrm{sign}(\tau)=1\end{subarray}}\prod_{j=1}^{n-k}b_{j\tau(j)}=\sum_{\begin{subarray}{c}\tau\in
D(n-k)\\\
\mathrm{sign}(\tau)=-1\end{subarray}}\prod_{j=1}^{n-k}b_{j\tau(j)}=0,$

  2. (2)

If $n-k$ is a even number, then

$\sum_{\begin{subarray}{c}\tau\in D(n-k)\\\
\mathrm{sign}(\tau)=(-1)^{\frac{n-k}{2}+1}\end{subarray}}\prod_{j=1}^{n-k}b_{j\tau(j)}=0.$

###### Proof.

First, let’s check this theorem when $\ell=n-k=1,2,3.$ If $\ell=1$, then $B=0$
and (1) holds; if $\ell=2$, then $B$ looks like

$\begin{pmatrix}0&\frac{1}{1-x^{-1}}\\\ \frac{1}{1-x}&0\end{pmatrix},x\neq 0.$

Since there is only a odd permutation in $D(2)$, thus (2) holds; if $\ell=3,$
then $B$ looks like

$xy\begin{pmatrix}0&\frac{1}{x-1}&\frac{1}{y-1}\\\
\frac{1}{1-x}&0&\frac{1}{y-x}\\\ \frac{1}{1-y}&\frac{1}{x-y}&0\\\
\end{pmatrix},x,y\neq 1\ \text{and}\ x\neq y.$

Since there are only two odd permutations $(1\ 2\ 3)$, $(1\ 3\ 2)$ in $D(3)$
and

$\frac{1}{x-1}\cdot\frac{1}{y-x}\cdot\frac{1}{1-y}+\frac{1}{y-1}\cdot\frac{1}{1-x}\cdot\frac{1}{x-y}=0.$

This implies that (1) holds

If $\ell\geqslant 4$ is odd, then

$\displaystyle\sum_{\begin{subarray}{c}\tau\in D(\ell)\\\
\mathrm{sign}(\tau)=1\end{subarray}}\prod_{j=1}^{\ell}b_{j\tau(j)}=x_{1}\cdots
x_{\ell}\sum_{\begin{subarray}{c}\tau\in D(\ell)\\\
\mathrm{sign}(\tau)=1\end{subarray}}\prod_{j=1}^{\ell}\frac{1}{x_{\tau(j)}-x_{j}},$

where $x_{1},\cdots,x_{\ell}$ are different from each other. For
$Y=\\{a_{1},\cdots,a_{r}\\}\subset X=\\{x_{1},\cdots,x_{\ell}\\}$, define

$f(Y)=\sum_{\tau\in L(r)}\prod_{j=1}^{r}\frac{1}{a_{\tau(r)}-a_{j}}.$

By lemma 3.1 we know that $f(X)=0$ when $|X|>2$. Note that

$\sum_{\begin{subarray}{c}\tau\in D(\ell)\\\
\mathrm{sign}(\tau)=1\end{subarray}}\prod_{j=1}^{\ell}\frac{1}{x_{\tau(j)}-x_{j}}=\sum_{\begin{subarray}{c}X=\bigsqcup_{i=1}^{s}X_{i}\\\
|X_{i}|\geqslant 2\\\ s\
\text{odd}\end{subarray}}\left(\prod_{j=1}^{s}f(X_{i})\right).$ (3.1)

Since $\ell$ is odd, every disjoint union $X=\bigsqcup_{i=1}^{s}X_{i}$ in the
RHS of (3.1) has an $X_{i}$ with $|X_{i}|>2$, thus

$\sum_{\begin{subarray}{c}\tau\in D(\ell)\\\
\mathrm{sign}(\tau)=1\end{subarray}}\prod_{j=1}^{\ell}\frac{1}{x_{\tau(j)}-x_{j}}=0.$

Similarly, it can be shown that $\displaystyle\sum_{\begin{subarray}{c}\tau\in
D(\ell)\\\
\mathrm{sign}(\tau)=-1\end{subarray}}\prod_{j=1}^{\ell}\frac{1}{x_{\tau(j)}-x_{j}}=0.$

If $\ell$ is even, one can prove (2) in the same way by using lemma 3.1. This
complete the proof of this theorem. ∎

Proof of equation (1.1) and equation (1.2):

In the notations of Theorem 3.2, let $x=\zeta$. If $k=0$ and $n$ is a even
number, then we have

$\sum_{\begin{subarray}{c}\tau\in D(n)\\\
\mathrm{sign}(\tau)=(-1)^{\frac{n}{2}+1}\end{subarray}}\prod_{j=1}^{n}a_{j\tau(j)}=0,$

which implies

$\sum_{\tau\in
D(n)}\prod_{j=1}^{n}a_{j\tau(j)}=(-1)^{\frac{n}{2}}\sum_{\tau\in
D(n)}\mathrm{sign}(\tau)\prod_{j=1}^{n}a_{j\tau(j)}=\dfrac{\left(\left(n-1\right)!!\right)^{2}}{2^{n}}.$

If $k=1,s_{k}=1$, and $n$ is an odd number, then we have

$\sum_{\begin{subarray}{c}\tau\in D(n-1)\\\
\mathrm{sign}(\tau)=(-1)^{\frac{n+1}{2}}\end{subarray}}\prod_{j=1}^{n-1}a_{j\tau(j)}=0,$

which implies that

$\sum_{\tau\in
D(n-1)}\prod_{j=1}^{n-1}a_{j\tau(j)}=(-1)^{\frac{n-1}{2}}\sum_{\tau\in
D(n)}\mathrm{sign}(\tau)\prod_{j=1}^{n+1}a_{j\tau(j)}=\dfrac{1}{n}\left(\dfrac{n-1}{2}!\right)^{2}.$

Thus we complete the proof of (1.1) and (1.2).∎

Acknowledgements: The author is deeply grateful to his colleagues Zhi-Wei Sun
and Keqin Liu for very helpful discussions.

## References

  * [1] F. Calogero and A.M. Perelomov, Some Diophantine relations involving circular functions of rational angles, Lin. Alg. Appl. 25, 91-94 (1979).
  * [2] K. Liu, On an identity from Sun’s conjectures, arxiv:submit/4338306.
  * [3] Z.-W. Sun, On some determinants with Legendre symbol entries, Finite Fields Appl. 56 (2019), 285-307.
  * [4] Z.-W. Sun, On some determinants involving the tangent function, preprint, arXiv:1901.04837.
  * [5] Z.-W. Sun, Permuations of $\\{1,\ldots,n\\}$ and related topics, J. Algebraic Combin., in press. See also arXiv:1811.10503.
  * [6] Z.-W. Sun, Arithmetic properties of some permanants, preprint, arxiv:2108.07723.
  * [7] H. Wang and Z.-W. Sun, Proof of a conjecture involving derangements and roots of unity, preprint.
  * [8] P. B. Denton, S. J. Parke, T. Tao, and X. Zhang, Eigenvectors from eigenvalues: a survey of a basic identity in linear algebra, Bulletin of the American Mathematical Society 59 (No. 1), 31-58.