Quantum Eigenvalue Estimation for Irreducible Non-negative Matrices
Anmer Daskin∗

arXiv:1505.02984v5 [quant-ph] 19 Nov 2015

Department of Computer Engineering, Istanbul Medeniyet University, Uskudar, Istanbul, Turkey
Quantum phase estimation algorithm has been successfully adapted as a sub frame of many other
algorithms applied to a wide variety of applications in different fields. However, the requirement of
a good approximate eigenvector given as an input to the algorithm hinders the application of the
algorithm to the problems where we do not have any prior knowledge about the eigenvector.
In this paper, we show that the principal eigenvalue of an irreducible non-negative operator can
be determined by using an equal superposition initial state in the phase estimation algorithm. This
removes the necessity of the existence of an initial good approximate eigenvector. Moreover, we
show that the success probability of the algorithm is related to the closeness of the operator to a
stochastic matrix. Therefore, we draw an estimate for the success probability by using the variance
of the column sums of the operator. This provides a priori information which can be used to know
the success probability of the algorithm beforehand for the non-negative matrices and apply the
algorithm only in cases when the estimated probability reasonably high. Finally, we discuss the
possible applications and show the results for random symmetric matrices and 3-local Hamiltonians
with non-negative off-diagonal elements.

I.

INTRODUCTION

In the recent decades, there have been many quantum algorithms proposed for the problems inefficient to solve on
classical computers. The quantum phase estimation algorithm [1, 2] is a leading illustration of such algorithms. It is
used for quantum simulations to estimate the eigenenergy corresponding to a given approximate eigenvector of the
unitary evolution operator of a quantum Hamiltonian. Furthermore, with different settings, it has been adapted as
a sub frame of many quantum algorithms applied to wide variety of applications in different fields (see the review
article ref.[3] and the references therein). However, the success of the algorithm partly depends on the pre-existing
approximate eigenvector. Thus, it may fail to output the right eigenvalue in the cases where the approximation to
the corresponding eigenvector is not good enough. This hinders the uses of the algorithm when there is not enough
prior knowledge to procure a good approximation to the eigenvector.
A symmetric matrix is called non-negative if all of its elements are greater than or equal to zero. It is irreducible
if the corresponding graph is strongly connected or if it is not reducible to a block-diagonal form by any column-row
permutation. Due to Perron-Frobenius theorem [4], an irreducible non-negative matrix have a positive eigenvector
(all elements are positive) with an associated positive principal eigenvalue whose magnitude is greater than the rest
of the eigenvalues. These matrices have been studied extensively [5, 6], the distribution of the coefficients of their
principal eigenvectors have been related to the matrix elements [7], and the sum of the coefficients have been shown
to be related to the number of walks in molecular graphs [8].
The phase estimation algorithm (PEA) mainly uses two quantum registers:|reg1i initially set to zero state and |reg2 i
holding an initial approximate eigenvector. In this article, we show that for irreducible non-negative matrices, one can
obtain the principal eigenvalue by using an equal superposition input state in PEA: i.e., instead of an approximate
eigenvector, initial value of |reg2 i is set to |0i and then put into equal superposition state by applying Hadamard
gates. In the output of the algorithm, this generates each one of the eigenvalues with the probability determined
by the normalized sum of the coefficients of the associated eigenvector. In addition, because all eigenvectors but the
principal one include positive-negative elements, we show that in most of the random cases the probability to see
principal eigenvalue in the output is much larger than the others. In also some cases, we show that first applying
Hadamard gates to the second register in the output (or measuring second register in the Hadamard basis) and then
measuring the first register when the measurement outcome of the second register is in |0 . . . 0i state increase the
success probability further in the output.
Eigenvector-coefficients of stochastic matrices sum to zero for all but the principal eigenvector. Therefore, for these
matrices one can generate the eigenvalue in the phase estimation algorithm with probability one [9]. Since any given
symmetric irreducible non-negative matrix can be converted into a stochastic matrix by a diagonal scaling matrix;
using the closeness of the matrix to a stochastic matrix, we also show that the success probability of the algorithm can

∗

email:anmerdaskin@yahoo.com

be predicted beforehand. We finally compare the estimated success probabilities with the computed ones for random
symmetric matrices and 3-local Hamiltonians of different dimensions with non-negative off-diagonal elements.
In the following sections, we shall first discuss PEA in detailed steps, then draw the estimates for the success
probability and finally show the possible applications of the algorithm.
II.

QUANTUM PHASE ESTIMATION ALGORITHM

The phase estimation algorithm (PEA) in general sense finds the value of φj for a given approximate eigenvector
|µj i in the eigenvalue equation U |µj i = eiφj |µj i. The algorithm mainly uses two quantum registers: Viz. |reg1i
initially set to zero state and |reg2 i holding an approximate eigenstate of a unitary matrix U . The first operation in
j
the algorithm puts |reg1i into the equal superposition state. In this setting, a sequence of operators, U 2 , controlled
by the jth qubit of |reg1i are applied to |reg2 i. Here, j = 1 . . . m and m is the number of qubits in |reg1i and also
determines the precision of the output. These sequential operations generate the quantum Fourier transform (QF T )
of the phase on |reg1i. Therefore, the application of the inverse quantum Fourier transform (QF T † ) turns the value
of |reg1i into the binary value of the phase. Consequently, one measures |reg1i to obtain the phase. Here, if the
unitary operator U is the time evolution operator of a quantum Hamiltonian H, i.e. U = eiH ; then one also obtains
the eigenenergy of that Hamiltonian.
For a symmetric irreducible non-negative matrix H of order 2n with ordered eigenvalues as φ1 ≥ · · · ≥ φ2n and
associated eigenvectors |µ1 i . . . |µ2n i, it is known that |µ1 i is the only eigenvector with all positive coefficients. Assume
j
unitary operator U = eiH and its powers U 2 with j = 0 . . . m are readily available for PEA. To estimate the value of
φ1 on the first register, in our setting, we simply modify the conventional phase estimation algorithm in the following
way:
• Instead of an approximate eigenvector, initial value of |reg2 i is set to |0i and then put into equal superposition
state by applying Hadamard gates.
• (optional) In the output, we change the basis of the second register by applying the Hadamard gates. Then,
if the measurement of the second register is equal to |0i state, then the principal eigenvalue is estimated on the
first register. This modification increases the success probability of the algorithm further when the probability
of measuring the principal eigenvalue is higher then the other eigenvalues.
In the following subsection, we shall describe the phase estimation algorithm in steps, also shown in Fig.1, to show
how the above modifications have an impact on the algorithm.
A.

Steps of the Algorithm

Here, the algorithm is described in details by showing how the quantum state changes in each step:
• The system is prepared in a way so that it includes two register: viz., |reg1 i and |reg2 i with m and n number
of qubits, respectively.
• Both registers are initialized into zero state: |ψ0 i=|reg1 i|reg2 i=|0i|0i.
• Then, the Hadamard operators are applied to both registers:
|ψ1 i = (H

⊗m

⊗H

⊗n

1

) |0i |0i = √
2n+m

2m+n
X−1
x=0

|xi

(1)
j

• As in the customary phase estimation algorithm; the unitary evolution operators U 2 controlled by the jth qubit
of |reg1 i are applied to |reg2 i, and finally the inverse of the quantum Fourier transform (QF T † ) is applied to
|reg1i.

At this point of the algorithm, we have a quantum state in which |reg1 i and |reg2 i hold the superposition of
the eigenvalues and the associated eigenvectors, respectively:
|ψ2 i =

N
−1
X
j=0

2

αj |φj i |µj i .

(2)

Here, the amplitude αj is related to the angle between the equal superposition state and the eigenvector |µj i and
P2n −1
can be described as the normalized sum of the coefficients of the eigenvectors since |reg2 i was √12n x=0 |xi at
the beginning:
n

2
√ X
αj = 1/ 2n
µji

(3)

i=1

p
If |reg2 i is written in the Hadamard basis, |+ · · · + +i+|+ · · · + −i+· · ·+|− · · · − −i, where |±i = 1/ (2)(|0i±
|1i); then the following quantum state is obtained:
n

|ψ2 i =

2
X
j=1

β1j |φj i |+ · · · +i
(4)

n

+ ···+

2
X
j=1

β2n j |φj i |− · · · −i

where βij s are new coefficients. Now, since there is only one eigenvector, viz. |µ1 i, with all positive real elements,
in the Hadamard basis, we can expect this positive eigenvector to be the closest state to |+ + · · · + +i.
• To revert the above state back to the standard basis, the Hadamard operator H ⊗n is again applied to |reg2 i:
|ψ3 i = (I ⊗ H ⊗n ) |ψ2 i

n

n

=

2
X
j=1

β1j |φj i |0 . . . 0i + · · · +

2
X
j=1

β2n j |φj i |1 . . . 1i

(5)

• |reg2 i is measured in the standard basis. As a result, for |reg2 i= |0 . . . 00i, the system collapses to the state
where the phase associated with the positive eigenvector is expected to be highly dominant in the first register.
• In the final step, the first register is measured to obtain the phase φ1 and get the eigenvalue of H. These steps
also drawn in Fig.1.

FIG. 1. Circuit design for the modified phase estimation algorithm. In the end, note that the application of the Hadamard
gates to the second register and the measurement on this register is optional.

B.

The Success Probability

It is easy to see that the success probability of the algorithm without the optional Hadamard gates on |reg2 i in
the output is determined from the normalized sum of the coefficients of the eigenvectors given in Eq.(3). Therefore,
the success probability for the dominant eigenvalue is α21 (note that the coefficients of |µ1 i are all positive, hence,
α21 = |α1 |2 .).
With the Hadamard gates in the output, we first find the probability to measure |reg2 i in |0i state and then find
the probability to see the dominant eigenvalue in |reg1 i knowing |reg2 i = |0i: As obvious in Eq.(5), the success
2
P n
2
. In addition, after measuring |0i on |reg2 i, the probability to measure
probability of getting |0i is
β
j=1 1j
3

|φ1 i on |reg1 i is |β11 |2 . Since we apply an equal superposition as an input, the component of an eigenvector in this
direction determines the probability to get the corresponding eigenvalue on |reg1 i. More formally, after the application
√ P2n
PN −1
of QF T † , we get |ψ2 i = j=0 αj |φj i |µj i, where αj = 1/ 2n i=1 µji . If |reg2 i of the state in Eq.(4) is measured in
the Hadamard basis, the probability to see |+ · · · +i is the sum of the amount of the components of the eigenvectors
in the direction of the initial vector:
n

Preg2 =

2
X
j=1

|αj |4 ≥

1
2n

(6)

Pn
This is also equal to 2j=1 β1j , i.e. the probability of measuring |reg2 i in |0 . . . 0i state at the end. Preg2 takes the
√
smallest value only when all |αj |s are equal to 1/ 2n . Moreover, when |reg2 i=|0 . . . 0i, the probability to see |φ1 i on
|reg1 i is:
Preg1 = |β11 |2 =

α41
.
Preg2

(7)

Although we have drawn the equalities for probabilities, without knowing the eigenvectors, it is not possible to
compute exact |αj |s and so Preg1 and Preg2 . Therefore,
√ we shall try to estimate them: Since α1 is the normalized sum
of the vector elements, it is easy to see that α1 ≥ 1/ 2n , where the equality occurs only when the magnitude of the
eigenvectors are the columns of the identity matrix. A similar observation is also made in ref.[9] where the principal
eigenvector is found for a given eigenvalue equal to 1.
Here, we will relate the matrix H to a stochastic matrix in order to develop some intuition for the estimation of the
probabilities. It is known that stochastic matrices have only one eigenvector with coefficients summing to a nonzero
value. A matrix can be made stochastic by a column wise scaling HD, where D is a diagonal matrix whose diagonal
elements are the inverses of the sums of the columns of H. The closeness of the matrix H to HD may provide a
prediction to estimate the success behavior of the phase estimation algorithm with an initial superposition state. The
closeness of H to HD can be defined by ||H − HD|| where || · || is any matrix norm. This also defines a perturbation
error. If we normalize this error term by ||H||, we get the relative error:
ǫ1 =

||H − HD||
≤ ||I − D||,
||H||

(8)

where I is an identity matrix. We can also look at the relative error in the inverses:
ǫ1 =

||H−1 − D−1 H−1 ||
≤ ||I − D−1 ||.
||H−1 ||

(9)

Since the matrix D also changes the eigenvector elements; when the variance of the diagonal elements of D is small,
we can expect the eigenvectors of H to show the similar behaviors as of HD and have one eigenvector whose elements
sum to a much greater value than the others. When H = HD or D = I, left and right principal eigenvectors of
symmetric H are the same and have the elements equal to √1N in the normalized space. Therefore, the variances σ1
and σ2 of the diagonal elements of D−1 and D, the column sums of H and their inverses, can somehow describe how
much the elements of the principal eigenvector deviate from √1N which is also the expectation value for randomly
generated normalized uniform N numbers. Using these intuitions, we define the estimate Preg2 as:
Λ1 + Λ2
Pereg2 =
2

(10)

with
Λ1 =

1
N
1
N

− σ1
, Λ2 =
+ σ1

1
N
1
N

− σ2
.
+ σ2

(11)

In addition, α21 = 1 for stochastic matrices and so for H when H = HD. Since the any deviation of the elements of
the principal eigenvector from √1N would change the equality α21 = 1, we multiply variance by N to get total deviation
from 1. Therefore, we define the estimate α21 as:
α
e21 =

(1 − N σ1 ) + (1 − N σ2 )
.
2
4

(12)

Note that α
e21 and Pereg2 are only predictions and may fail to provide good estimates for α21 and Preg2 , respectively,
for particularly some structured matrices and matrices whose eigenvectors are close to the columns of an identity
matrix. For instance, consider the following matrix:


21.8214
0
0.6118 0.4983
14.2944 0.4983 0.6118
 0
.
 0.6118 0.4983 12.1626
0 
0.4983 0.6118
0
5.4111

(13)

The eigenvalues of this matrix are 5.3537, 12.0193, 14.4411, and 21.8753, respectively, associated with the following
eigenvectors:


0.0304613
 0.0686662
−0.0077623
−0.9971443

0.0597207
0.2074209
−0.9761393
0.0237068

0.0215934
−0.9758165
−0.2076080
−0.0649217


0.9975166
0.0066086
.
0.0631720
0.0304360

(14)

As computed from the above, the magnitudes of the sums of the components of the eigenvectors are very close to
each other: | − 0.90578|, | − 0.68529|, | − 1.22675|, and |1.09773|, respectively. In addition, the second eigenvector, not
associated to the principal eigenvalue, has the largest sum.
Also note that one can also predict α21 from the distribution of the matrix elements. For instance [7, 10],
maxij

µ1i
maxij hij
≤
,
µ1j
minij hij

(15)

where hij s are the nonzero matrix elements of H. However, these bounds are generally not tight enough to give good
estimates in many cases.
In Sec.IV, we shall show the comparisons of these estimates with the computed probabilities for random matrices.
III.

POSSIBLE APPLICATIONS

In the previous sections, we have showed the positive eigenpair of an irreducible symmetric non-negative matrix can
be obtained by using an initial equal superposition state when the estimated probabilities high. This eliminates the
necessity of knowing an initial approximate eigenvector in the applications of the phase estimation algorithm to the
problems involving non-negative matrices. Irreducible non-negative matrices are known to have positive eigenpairs
due to Perron-Frobenious theorem. There is also more general but the same statements indicated in this theorem for
compact operators [11]. One may also apply permutation matrices or splitting techniques to convert the interested
matrix to a non-negative one or to increase the expected probability outcome when the probability is low than an
acceptable value.
Since a wide variety of problems in science can be represented by non-negative matrices, the phase estimation
algorithm can be applied without an initial approximate eigenvector to these problems. In the following subsections,
we shall consider two important classes among these problems: viz., the k-local Hamiltonian problem and the stochastic
processes.
A.

Applications to Local Hamiltonian Problems

One of the applications of the non-negative matrices can be found in quantum mechanics where such matrices called
stoquastic matrices [12, 13]: i.e., matrices with non-negative off-diagonal elements. If a Hamiltonian operator H, a
Hermitian operator acting on (C 2 )⊗n , can be expressed as a sum of local Hamiltonians representing
P the interactions
between at most k qubits, then it is called k-local n-qubit Hamiltonian. More formally, H = s Hs . Finding the
lowest eigenvalue of H defines an eigenvalue problem which is so called local Hamiltonian problem. This eigenvalue
problem can be also described as a decision problem to decide whether the ground state energy of H is at most a
or at least b for given constants a and b such that a ≤ b and b − a ≤ 1/poly(n). This problem has been shown to
be QM A-complete for k ≥ 2 [14] (for k=2, it is QMA-complete only when both negative and positive signs in the
Hamiltonian are allowed [15]).
For the local Hamiltonians with positive off-diagonal matrix elements in the standard basis, the matrix elements
of the corresponding Gibbs density matrix, ρ = e−βH /T r(−βH), are all positive for any β ≥ 0. Due to the theorem
5

given above, The ground state energy of H with positive off-diagonal elements can be associated with the eigenvector
whose coefficients are all positive. Because one can associate a probability distribution to the ground state, the nature
of these Hamiltonians are considered similar to stochastic processes. Thus, the Hamiltonians of these types are called
stoquastic [12].(note that the Hamiltonian is not a stochastic matrix whose rows and columns are sum to one and the
principal eigenvalue is one associated with an eigenvector of all ones.)
B.

Applications to Stochastic Processes

The phase estimation algorithm can find the principal eigenpair of a stochastic matrix with probability one [9] since
the sum of the coefficients is one for the principal eigenvector and zero for the rest of the eigenvectors. However,
the stochastic processes are generally defined by non-Hermitian matrices, which are also widely seen in quantum
physics: The Hamiltonian of a closed system in equilibrium must be Hermitian to have real energy values. However,
nonequilibrium processes and open systems connected to reservoirs can be only described by non-Hermitian models
[16].
1.

Simulation of Non-Hermitian Operators

Any non-Hermitian matrix H can be decomposed into a Hermitian and a skew-Hermitian parts as:
H=H +S =

1
1
(H + H† ) + (H − H† )
2
2

(16)

where H† describes the conjugate transpose of H. Here, H = 21 (H + H† ) and S = 21 (H − H† ) define the nearest Hermitian and skew Hermitian matrices to H, respectively [17].Therefore, the matrix H can be used as an approximation
to H in the simulation.
A matrix H is called normal if H† H − HH† = 0. When H is a normal, it is easy to see that H and S also commute:
[H, S] = HS − HS = 0. Moreover, because of the symmetry, all the eigenvalues of H are real; and because of the
skew-symmetry, all the eigenvalues of S involve only imaginary parts. Since HH† = H† H, the eigenvectors of H and
S are the same. In addition, the imaginary part of the eigenvalues of H are equal to the eigenvalues of S and the
real parts are to the eigenvalues of H. Hence, using U1 = eiH and U2 = eS in the simulation, one can simulate the
non-Hermitian operator H on quantum computers by using separate two registers to obtain the imaginary and the
real parts of the eigenvalue. However, in the stochastic cases, if H is normal, it must be also doubly stochastic: i.e., its
left and right principal eigenvectors are already known to be a vector of all ones with the eigenvalue one. Therefore,
instead of an approximate normal matrix, we can use the closest Hermitian matrix defined above as H = 12 (H + H† ).
In the non-stochastic cases, one can also try to instead find the closest normal matrix by following the Jacobi
algorithm which attempts to diagonalize the matrix by using rotation matrices (Givens rotations). and converges to
a matrix so called ∆H. The best closest matrix in the Frobenius norm is then defined by the sum of the diagonal
elements of ∆H and the rotation matrices used in the algorithm [18, 19]. One can also procure a quantum circuit
to implement this closest matrix by mapping rotation matrices to quantum gates as done in [20]. However, since
this algorithm is based on the eigenvalue decomposition, the eigenvalues (the diagonal elements of ∆H, and the
eigenvectors (the combination of the rotation matrices ) are already generated for the found normal matrix.
IV.

NUMERICAL RESULTS

e21 , Pereg1 , Pereg2 , respectively, with the computed
In this section, we compare the estimated α21 , Preg1 , Preg2 : i.e. α
2
probabilities α1 , Preg1 , Preg2 for the randomly generated matrices described below.
A.

Random Symmetric Matrices with Non-Negative Off-Diagonal Elements

In MATLAB, we generate a random symmetric matrix with non-negative off-diagonal elements as follows: first
a random diagonal matrix and a strictly upper triangular sparse non-negative matrix are generated. Then, these
matrices are combined to have a symmetric matrix with non-negative off-diagonal elements:
%MATLAB code to generate a random symmetric matrix
% with nonnegative off-diagonal elements
6

R
%
L
H

= triu(sprand(N,N,0.5),1); % N is the dimension,
and 0.5 is the sparsity of the matrix.
= randn(N,1);
= R + R’ + diag(L);

Fig.2 shows the comparisons of the estimated and computed probabilities for matrices of different dimensions generated
by following the above code. As shown in the figure, the estimated and computed probabilities are almost the same
and very high (close to 1). Furthermore, the success probability, Preg1 , with the Hadamard gates applied to |reg2 i in
the end is almost one and higher than the probability, α21 , without the Hadamard gates.
B.

Random 3-Local Hamiltonians

As our first example, the following Hamiltonian is employed:
X
HXZ =
Kijk Xi Xj Xk + Jijk Zi Zj Zk

(17)

i,j,k

with 0 ≤ Kijk ≤ 1 and −1 ≤ Jijk ≤ 1. Here, X and Z are the Pauli spin operators. Choosing Kijk and Jijk randomly,
50 numbers of random matrices are generated for each 9, 10, 11, 12, and 13 qubit systems and show the estimated
and computed probabilities in Fig.3.
In the second example, 1-body and 2-body interaction terms are also included separately in the Hamiltonian:
X
X
K1ij Xi Xj + J1ij Zi Zj
di Xi + hi Zi +
HXZ =
i,j

i

+

X

K2ijk Xi Xj Xk + J2ijk Zi Zj Zk

(18)

i,j,k

with 0 ≤ di , K1ijk , K2ijk ≤ 1 and −1 ≤ hi , J1ijk , J2ijk ≤ 1. As done for Eq.(17, choosing the coefficients randomly,
Hamiltonians of different sizes are generated. The results are shown in Fig.4. As seen in the figure, the estimation is
not as good as the cases in Fig.2 and Fig.3. This is because the ratio defined in Eq.(15) is generally higher for the
matrices generated from Eq.(18) in comparison to the ones from Eq.(17).
V.

CONCLUSION

In this paper, we have showed the positive eigenpair of an irreducible symmetric non-negative matrix can be obtained
when the estimated probabilities high. This eliminates the necessity of knowing an initial approximate eigenvector
in the applications of the phase estimation algorithm to a wide variety of problems involving non-negative matrices.
Here, we also have discussed how to apply it to local Hamiltonian problems and stochastic processes and showed
the comparisons of the estimated success probabilities with the computed ones for random symmetric matrices and
3-local Hamiltonians.

[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]

D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 83, 5162 (1999).
A. Kitaev, Electronic Colloquium on Computational Complexity (ECCC) 3 (1996).
I. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153 (2014).
C. Meyer, Matrix analysis and applied linear algebra book and solutions manual, Vol. 2 (Society for Industrial and Applied
Mathematics, 2000).
R. B. Bapat and T. E. Raghavan, Nonnegative matrices and applications, Vol. 64 (Cambridge University Press, 1997).
A. Berman and R. J. Plemmons, The Mathematical Sciences, Classics in Applied Mathematics 9 (1979).
H. Minc, SIAM Journal on Numerical Analysis 7, 424 (1970).
I. Gutman, C. Rucker, and G. Rucker, Journal of Chemical Information and Computer Sciences 41, 739 (2001).
A. Daskin, A. Grama, and S. Kais, Quantum Information Processing 13, 2653 (2014).
S. M. Cioaba and D. A. Gregory, Electron. J. Linear Algebra 16, 366 (2007).
Y. Du, Order structure and topological methods in nonlinear partial differential equations, Vol. 2 (World Scientific, 2006)
Chap. 1: Krein-Rutman Theorem and the Principal Eigenvalue.
S. Bravyi, D. P. Divincenzo, R. Oliveira, and B. M. Terhal, Quantum Info. Comput. 8, 361 (2008).

7

[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]

S. Bravyi and B. Terhal, SIAM Journal on Computing 39, 1462 (2010), http://dx.doi.org/10.1137/08072689X.
J. Kempe, A. Kitaev, and O. Regev, SIAM Journal on Computing 35, 1070 (2006).
S. P. Jordan, D. Gosset, and P. J. Love, Phys. Rev. A 81, 032331 (2010).
K. B. Efetov, Physical review letters 79, 491 (1997).
J. B. Keller, Mathematics Magazine 48, pp. 192 (1975).
N. J. Higham, in Applications of Matrix Theory (University Press, 1989) pp. 1–27.
A. Ruhe, BIT Numerical Mathematics 27, 585 (1987).
J. J. Vartiainen, M. Möttönen, and M. M. Salomaa, Physical review letters 92, 177902 (2004).

0.996
Preg2
Pereg2

Preg2
Pereg2

0.992

Probability

Probability

0.993

0.991
0.99

0.9955

0.995

0.9945

0.989

1.001
Preg1
Pereg1
α 12
α
e 12

1

Probability

Probability

1.002

0.998
0.996
0.994
0

5

10

15

20

25
Matrix

30

35

40

45

0.999
0.998
0.997
0

50

Preg1
Pereg1
α 12
α
e 12

1

5

(a)512 by 512 Matrices

25
Matrix

30

35

40

45

50

Preg2
Pereg2

0.9989
Probability

Probability

20

0.9989
Preg2
Pereg2

0.9977
0.9976
0.9975

0.9988
0.9988
0.9988

0.9974

0.9988

1

1
Preg1
Pereg1
α 12
α
e 12

0.9995

Probability

Probability

15

(b)1024 by 1024 Matrices

0.9978

0.999

0.9985
0

10

5

10

15

20

25
Matrix

30

35

40

45

0.9996
0.9994
0.9992
0

50

(c)2048 by 2048 Matrices

Preg1
Pereg1
α 12
α
e 12

0.9998

5

10

15

20

α4
1
Preg2

30

35

40

45

50

(d)4096 by 4096 Matrices

FIG. 2. Representations of the comparisons of the estimated probabilities α
e21 , Pereg1 =
probabilities α21 , Preg1 =

25
Matrix

α
e4
1
ereg
P
2

, Pereg2 with the computed

, and Preg2 for random non-negative symmetric matrices of different dimensions.

8

0.96

0.95
Preg2
Pereg2
Probability

Probability

0.9
0.85
0.8
0.75

0.9

1

1
Preg1
Pereg1
α 12
α
e 12

0.95

Probability

Probability

0.92

0.88

0.7

0.9
0.85
0.8
0

Preg2
Pereg2

0.94

5

10

15

20

25
Matrix

30

35

40

45

0.96

0.94
0

50

Preg1
Pereg1
α 12
α
e 12

0.98

(a)512 by 512 3-local Hamiltonians

5

10

15

20

25
Matrix

30

35

40

45

(b)1024 by 1024 3-local Hamiltonians

0.98

0.98
Preg2
Pereg2
Probability

Probability

Preg2
Pereg2
0.96

0.94

0.92

0.97

0.96

0.95

1.01

1

0.99

Probability

Preg1
Pereg1
α 12
α
e 12

1
Probability

50

0.98

Preg1
Pereg1
α 12
α
e 12

0.99

0.98

0.97
0.96
0

5

10

15

20

25
Matrix

30

35

40

45

0.97
0

50

(c)2048 by 2048 3-local Hamiltonians

5

10

15

20

probabilities α21 , Preg1 =

30

35

40

45

50

(d)4096 by 4096 3-local Hamiltonians

FIG. 3. Representations of the comparisons of the estimated probabilities α
e21 , Pereg1 =
α4
1
Preg2

25
Matrix

α
e4
1
ereg
P
2

, Pereg2 with the computed

, and Preg2 for random 3-local Hamiltonians of different dimensions described in Eq.(17).

9

1

1

Preg2
Pereg
Probability

2

0.9
0.85

Probability

0.95

0.9

Preg2
Pereg2

0.8
0.75

0.85

1

1

Probability

Probability

0.95

0.95
Preg1
Pereg1
α 12
α
e 12

0.9

0.85
0

5

10

15

20

25
Matrix

30

35

40

45

0.96
0.94
0.92
0

50

Preg1
Pereg1
α 12
α
e 12

0.98

(a)512 by 512 3-local Hamiltonians

5

10

15

20

25
Matrix

30

35

40

45

50

(b)1024 by 1024 3-local Hamiltonians

0.98

1

Probability

Probability

0.96
0.94
0.92

Preg2
Pereg2

0.88

0.85

1

1
Preg1
Pereg1
α 12
α
e 12

0.98

Probability

Probability

0.9

Preg2
Pereg2

0.9

0.96

0.94
0

0.95

5

10

15

20

25
Matrix

30

35

40

45

0.98
0.96

0.92
0

50

(c)2048 by 2048 3-local Hamiltonians

Preg1
Pereg1
α 12
α
e 12

0.94

5

10

15

20

probabilities α21 , Preg1 =

30

35

40

45

50

(d)4096 by 4096 3-local Hamiltonians

FIG. 4. Representations of the comparisons of the estimated probabilities α
e21 , Pereg1 =
α4
1
Preg2

25
Matrix

α
e4
1
ereg
P
2

, Pereg2 with the computed

, and Preg2 for random 3-local Hamiltonians of different dimensions described in Eq.(18).

10

