# Computing Coherence Vectors and Correlation Matrices with Application to Quantum Discord Quantification

Authors: Jonas Maziero
Journal: Advances in Mathematical Physics
 (2016)

---

## Abstract

Coherence vectors and correlation matrices are important functions frequently used in physics. The numerical calculation of these functions directly from their definitions, which involves Kronecker products and matrix multiplications, may seem to be a reasonable option. Notwithstanding, as we demonstrate in this paper, some algebraic manipulations before programming can reduce considerably their computational complexity. Besides, we provide Fortran code to generate generalized Gell-Mann matrices and to compute the optimized and unoptimized versions of associated Bloch’s vectors and correlation matrix in the case of bipartite quantum systems. As a code test and application example, we consider the calculation of Hilbert-Schmidt quantum discords.


---

## Body

## 1. Introduction

Correlation functions are fundamental objects for statistical analysis and are thus ubiquitous in most kinds of scientific inquiries and their applications [1, 2]. In physics, correlation functions have an important role for research in areas such as quantum optics and open systems [3, 4], phase transitions and condensed matter physics [5, 6], and quantum field theory and nuclear and particle physics [7]. Another area in which correlation functions are omnipresent is quantum information science (QIS), an interdisciplinary field that extends the applicability of the classical theories of information, computation, and computational complexity [8–10].Investigations about the quantum correlations in physical systems have been one of the main catalyzers for developments in QIS [11–15]. There are several guises of quantum correlations, and quantum discord stands among the most promising quantum resources for fueling the quantum advantage [16–26]. When computing or witnessing quantum discord, or other kinds of correlation or quantumness quantifiers, we are frequently faced with the need for calculating coherence vectors and correlation matrices [27–41]. And it is the main aim of this paper to provide formulas for these functions which are amenable for more efficient numerical calculations when compared with the direct implementation of their definitions.In order to define coherence vectors and correlation matrices, let us consider a composite bipartite system with Hilbert spaceH
a
b
=
H
a
⊗
H
b. Hereafter, the corresponding dimensions are denoted by d
s
=
dim
⁡
H
s for s
=
a
b
,
a
,
b. In addition, let Γ
j
s, with(1)
T
r
Γ
j
s
=
0
,
T
r
Γ
j
s
Γ
k
s
=
2
δ
j
kbeing a basis for the special unitary group S
U
(
d
s
). Any density operator describing the state of the system H
a
b can be written in the local basis Γ
j
a
⊗
Γ
k
b as follows:(2)
ρ
=
1
d
a
b
I
a
⊗
I
b
+
∑
j
=
1
d
a
2
-
1
a
j
Γ
j
a
⊗
I
b
+
I
a
⊗
∑
k
=
1
d
b
2
-
1
b
k
Γ
k
b
+
∑
j
=
1
d
a
2
-
1
∑
k
=
1
d
b
2
-
1
c
j
,
k
Γ
j
a
⊗
Γ
k
b
,where(3)
j
=
1
,
…
,
d
a
2
-
1
,
k
=
1
,
…
,
d
b
2
-
1and I
s is the identity operator in H
s. One can readily verify that the components of the coherence (or Bloch’s) vectors a
=
(
a
1
,
…
,
a
d
a
2
-
1
) and b
=
(
b
1
,
…
,
b
d
b
2
-
1
) and of the correlation matrix C
=
(
c
j
,
k
) are given by(4)
a
j
=
2
-
1
d
a
T
r
Γ
j
a
⊗
I
b
ρ
,
b
k
=
2
-
1
d
b
T
r
I
a
⊗
Γ
k
b
ρ
,
c
j
,
k
=
2
-
2
d
a
b
T
r
Γ
j
a
⊗
Γ
k
b
ρ
.It is worth mentioning that the mean value ofany observable in H
s, for s
=
a
,
b
,
a
b, can be obtained using these quantities.Inhttps://github.com/jonasmaziero/LibForQ.git, we provide Fortran code to compute the coherence vectors, correlation matrices, and quantum discord quantifiers we deal with here. Besides these functions, there are other tools therein that may be of interest to the reader. The instructions on how to use the software are provided in the readme file. Related to the content of this section, the subroutine bloch_vector_gellmann_unopt(d
s, ρ
s, s) returns the coherence vectors a or b and the subroutine corrmat_gellmann_unopt(d
a, d
b, ρ, C) computes the correlation matrix C. Now, let us notice that if calculated directly from the equations above, for d
a
,
d
b
≫
1, the computational complexity (CC) to obtain the coherence vectors a and b or the correlation matrix C is(5)
C
C
a
=
C
C
b
=
C
C
C
≈
O
d
a
6
d
b
6
.The remainder of this paper is structured as follows. In Section2, we obtain formulas for a, b, and C which are amenable for more efficient numerical computations. In Section 3, we test these formulas by applying them in the calculation of Hilbert-Schmidt quantum discords. In Section 4, we make some final remarks about the usefulness and possible applications of the results reported here.


## 2. Computing Coherence Vectors and Correlation Matrices

The partial trace function [42] can be used in order to obtain the reduced states ρ
a
=
T
r
b
(
ρ
) and ρ
b
=
T
r
a
(
ρ
) and to write the components of Bloch’s vectors in the form(6)
a
j
=
2
-
1
d
a
T
r
Γ
j
a
ρ
a
,
b
k
=
2
-
1
d
b
T
r
Γ
k
b
ρ
b
.Thus, when computing the coherence vectors of the parties a and b, we shall solve a similar problem; so let us consider it separately. That is to say, we shall regard a generic density operator written as(7)
ρ
s
=
1
d
s
I
s
+
∑
j
=
1
d
2
-
1
s
j
Γ
j
s
,where s
j
=
2
-
1
d
s
T
r
(
ρ
s
Γ
j
s
).Now, and for the remainder of this paper, we assume that the matrix elements of regarded density operatorρ in the standard computational basis are given. We want to computeBloch’s vector [43]: s
=
(
s
1
,
…
,
s
d
s
2
-
1
). For the sake of doing that, a particular basis Γ
j
s must be chosen. Here we pick the generalized Gell-Mann matrices, which are presented as follows in three groups [44]:(8)
Γ
j
s
1
=
2
j
j
+
1
∑
k
=
1
j
+
1
-
j
δ
k
,
j
+
1
k
k
,
for
j
=
1
,
…
,
d
s
-
1
,
Γ
k
,
l
s
2
=
k
l
+
l
k
,
for
1
≤
k
<
l
≤
d
s
,
Γ
k
,
l
s
3
=
-
i
k
l
-
l
k
,
for
1
≤
k
<
l
≤
d
s
,which are named the diagonal, symmetric, and antisymmetric groups, respectively. The last two groups possess d
s
(
d
s
-
1
)
/
2 generators each. Any one of these matrices can be obtained by calling the subroutine gellmann(d
s, g, k, l, Γ
(
k
,
l
)
s
(
g
)). For the first group, g
=
1, we make j
=
k and, in this case, one can set l to any integer.It is straightforward seeing that, for the generators above, the corresponding components of Bloch’s vector can be expressed directly in terms of the matrix elements of the density operatorρ
s as follows:(9)
s
j
1
=
d
s
2
j
j
+
1
∑
k
=
1
j
+
1
-
j
δ
k
,
j
+
1
k
ρ
s
k
,
for
j
=
1
,
…
,
d
s
-
1
,
s
k
,
l
2
=
d
s
R
e
l
ρ
s
k
,
for
1
≤
k
<
l
≤
d
s
,
s
k
,
l
2
=
d
s
I
m
l
ρ
s
k
,
for
1
≤
k
<
l
≤
d
s
.These expressions were implemented in the Fortran subroutinebloch_vector_gellmann(d
s, ρ
s, s). With this subroutine and the partial trace function [42], we can compute the coherence vectors a and b.We observe that after these simple algebraic manipulations the computational complexity of Bloch’s vector turns out to be basically the CC for the partial trace function. Hence, from [42], we have that, for d
a
,
d
b
≫
1,(10)
C
C
a
≈
O
d
a
2
d
b
,
C
C
b
≈
O
d
a
d
b
2
.One detail we should keep in mind when making use of the codes linked to this paper is the convention we apply to the indexes of the components ofs. For the first group of generators, Γ
j
s
(
1
), naturally, j
=
1
,
…
,
d
s
-
1. We continue with the second group of generators, Γ
j
s
(
2
)
=
Γ
(
k
,
l
)
s
(
2
), by setting j
(
k
,
l
)
=
(
1,2
)
=
d
s
-
1
+
1
=
d
s, j
(
k
,
l
)
=
(
1,3
)
=
d
s
+
1
,
…
,
j
(
k
,
l
)
=
(
1
,
d
s
)
=
2
(
d
s
-
1
), j
(
k
,
l
)
=
(
2,3
)
=
2
(
d
s
-
1
)
+
1
,
…. The same convention is used for the third group of generators, Γ
j
s
(
3
)
=
Γ
(
k
,
l
)
s
(
3
), but here we begin with j
(
k
,
l
)
=
(
1,2
)
=
d
s
-
1
+
2
-
1
d
s
(
d
s
-
1
)
+
1
=
d
s
+
2
-
1
d
s
(
d
s
-
1
).Next we address the computation of thecorrelation matrix
C
=
(
c
j
,
k
), which is a d
a
2
-
1
×
(
d
b
2
-
1
) matrix that we write in the following form:(11)
C
=
C
1,1
C
1,2
C
1,3
C
2,1
C
2,2
C
2,3
C
3,1
C
3,2
C
3,3
,with the submatrices given as shown below. For convenience, we define the auxiliary variables:(12)
ι
≔
2
j
j
+
1
,
κ
≔
2
k
k
+
1
,
ς
≔
d
a
d
b
4
.The matrix elements ofC
(
1,1
), whose dimension is d
a
-
1
×
(
d
b
-
1
), correspond to the diagonal generators for a and diagonal generators for b:(13)
c
j
,
k
1,1
=
ς
T
r
Γ
j
a
1
⊗
Γ
k
b
1
ρ
=
ς
ι
κ
T
r
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
m
m
⊗
∑
p
=
1
k
+
1
-
k
δ
p
,
k
+
1
p
p
ρ
=
ς
ι
κ
∑
m
=
1
j
+
1
∑
p
=
1
k
+
1
-
j
δ
m
,
j
+
1
-
k
δ
p
,
k
+
1
T
r
m
m
⊗
p
p
ρ
=
ς
ι
κ
∑
m
=
1
j
+
1
∑
p
=
1
k
+
1
-
j
δ
m
,
j
+
1
-
k
δ
p
,
k
+
1
m
p
ρ
m
p
.The matrix elements ofC
(
1,2
), whose dimension is d
a
-
1
×
2
-
1
d
b
(
d
b
-
1
), correspond to the diagonal generators for a and symmetric generators for b:(14)
c
j
,
k
1,2
=
ς
T
r
Γ
j
a
1
⊗
Γ
k
b
2
ρ
=
ς
T
r
Γ
j
a
1
⊗
Γ
p
,
q
b
2
ρ
=
ς
T
r
ι
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
m
m
⊗
p
q
+
q
p
ρ
=
ς
ι
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
T
r
m
p
m
q
ρ
+
T
r
m
q
m
p
ρ
=
ς
ι
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
m
q
ρ
m
p
+
m
q
ρ
m
p
∗
=
2
ς
ι
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
R
e
m
q
ρ
m
p
.The matrix elements ofC
(
1,3
), whose dimension is d
a
-
1
×
2
-
1
d
b
(
d
b
-
1
), correspond to the diagonal generators for a and antisymmetric generators for b:(15)
c
j
,
k
1,3
=
ς
T
r
Γ
j
a
1
⊗
Γ
k
b
3
ρ
=
ς
T
r
Γ
j
a
1
⊗
Γ
p
,
q
b
3
ρ
=
-
i
ς
ι
T
r
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
m
m
⊗
p
q
-
q
p
ρ
=
-
i
ς
ι
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
T
r
m
p
m
q
ρ
-
T
r
m
q
m
p
ρ
=
-
i
ς
ι
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
m
q
ρ
m
p
-
m
q
ρ
m
p
∗
=
2
ς
ι
∑
m
=
1
j
+
1
-
j
δ
m
,
j
+
1
I
m
m
q
ρ
m
p
.The matrix elements ofC
(
2,1
), whose dimension is 2
-
1
d
a
d
a
-
1
×
(
d
b
-
1
), correspond to the symmetric generators for a and diagonal generators for b:(16)
c
j
,
k
2,1
=
ς
T
r
Γ
j
a
2
⊗
Γ
k
b
1
ρ
=
ς
T
r
Γ
m
,
n
a
2
⊗
Γ
k
b
1
ρ
=
ς
T
r
m
n
+
n
m
⊗
κ
∑
p
=
1
k
+
1
-
k
δ
p
,
k
+
1
p
p
ρ
=
ς
κ
∑
p
=
1
k
+
1
-
k
δ
p
,
k
+
1
T
r
m
p
n
p
ρ
+
T
r
n
p
m
p
ρ
=
ς
κ
∑
p
=
1
k
+
1
-
k
δ
p
,
k
+
1
n
p
ρ
m
p
+
n
p
ρ
m
p
∗
=
2
ς
κ
∑
p
=
1
k
+
1
-
k
δ
p
,
k
+
1
R
e
n
p
ρ
m
p
.The matrix elements ofC
(
2,2
), whose dimension is 2
-
1
d
a
d
a
-
1
×
2
-
1
d
b
(
d
b
-
1
), correspond to the symmetric generators for a and symmetric generators for b:(17)
c
j
,
k
2,2
=
ς
T
r
Γ
j
a
2
⊗
Γ
k
b
2
ρ
=
ς
T
r
Γ
m
,
n
a
2
⊗
Γ
p
,
q
b
2
ρ
=
ς
T
r
m
n
+
n
m
⊗
p
q
+
q
p
ρ
=
ς
n
q
ρ
m
p
+
m
q
ρ
n
p
+
n
p
ρ
m
q
+
m
p
ρ
n
q
=
2
ς
R
e
n
q
ρ
m
p
+
R
e
n
p
ρ
m
q
.The matrix elements ofC
(
2,3
), whose dimension is 2
-
1
d
a
d
a
-
1
×
2
-
1
d
b
(
d
b
-
1
), correspond to the symmetric generators for a and antisymmetric generators for b:(18)
c
j
,
k
2,3
=
ς
T
r
Γ
j
a
2
⊗
Γ
k
b
3
ρ
=
ς
T
r
Γ
m
,
n
a
2
⊗
Γ
p
,
q
b
3
ρ
=
ς
T
r
m
n
+
n
m
⊗
-
i
p
q
-
q
p
ρ
=
-
i
ς
n
q
ρ
m
p
+
m
q
ρ
n
p
-
n
p
ρ
m
q
-
m
p
ρ
n
q
=
2
ς
I
m
n
q
ρ
m
p
-
I
m
n
p
ρ
m
q
.The matrix elements ofC
(
3,1
), whose dimension is 2
-
1
d
a
d
a
-
1
×
(
d
b
-
1
), correspond to the antisymmetric generators for a and diagonal generators for b:(19)
c
j
,
k
3,1
=
ς
T
r
Γ
j
a
3
⊗
Γ
k
b
1
ρ
=
T
r
Γ
m
,
n
a
3
⊗
Γ
k
b
1
ρ
=
ς
T
r
-
i
m
n
-
n
m
⊗
κ
∑
p
=
1
k
+
1
-
k
δ
p
,
k
+
1
p
p
ρ
=
-
i
ς
κ
∑
p
=
1
k
+
1
-
k
δ
p
,
k
+
1
n
p
ρ
m
p
-
m
p
ρ
n
p
=
2
ς
κ
∑
p
=
1
k
+
1
-
k
δ
p
,
k
+
1
I
m
n
p
ρ
m
p
.The matrix elements ofC
(
3,2
), whose dimension is 2
-
1
d
a
d
a
-
1
×
2
-
1
d
b
(
d
b
-
1
), correspond to the antisymmetric generators for a and symmetric generators for b:(20)
c
j
,
k
3,2
=
ς
T
r
Γ
j
a
3
⊗
Γ
k
b
2
ρ
=
ς
T
r
Γ
m
,
n
a
3
⊗
Γ
p
,
q
b
2
ρ
=
ς
T
r
-
i
m
n
-
n
m
⊗
p
q
+
q
p
ρ
=
-
i
ς
n
q
ρ
m
p
-
m
p
ρ
n
q
+
n
p
ρ
m
q
-
m
q
ρ
n
p
=
2
ς
I
m
n
q
ρ
m
p
+
I
m
n
p
ρ
m
q
.The matrix elements ofC
(
3,3
), whose dimension is 2
-
1
d
a
d
a
-
1
×
2
-
1
d
b
(
d
b
-
1
), correspond to antisymmetric generators for a and antisymmetric generators for b:(21)
c
j
,
k
3,3
=
ς
T
r
Γ
j
a
3
⊗
Γ
k
b
3
ρ
=
ς
T
r
Γ
m
,
n
a
3
⊗
Γ
p
,
q
b
3
ρ
=
ς
T
r
-
i
m
n
-
n
m
⊗
-
i
p
q
-
q
p
ρ
=
-
ς
n
q
ρ
m
p
+
m
p
ρ
n
q
-
n
p
ρ
m
q
-
m
q
ρ
n
p
=
2
ς
R
e
n
p
ρ
m
q
-
R
e
n
q
ρ
m
p
.We remark that when implementing these expressions numerically, for the sake of mapping the local to the global computational basis, we utilize, for example,(22)
n
p
≡
n
-
1
d
b
+
p
.The subroutinecorrmat_gellmann(d
a, d
b, ρ, C) returns the correlation matrix C
=
(
c
j
,
k
), as written in (11), associated with the bipartite density operator ρ and computed using the Gell-Mann basis, as described in this section. The convention for the indexes of the matrix elements c
j
,
k is defined in the same way as for the coherence vectors. The computational complexity for C, computed via the optimized expressions obtained in this section, is, for d
a
,
d
b
≫
1,(23)
C
C
C
≈
O
d
a
2
d
b
2
.By generating some random density matrices [45], we checked that the expressions and the corresponding code for the unoptimized and optimized versions of a, b, and C agree. Additional tests shall be presented in the next section, where we calculate some quantum discord quantifiers.


## 3. Computing Hilbert-Schmidt Quantum Discords

The calculation of quantum discord functions (QD) usually involves hard optimization problems [46, 47]. In the last few years, a large amount of effort has been dedicated towards computing QD analytically, with success being obtained mostly for low-dimensional quantum systems [48–65]. Although not meeting all the required properties for abona fide QD quantifier [66], the Hilbert-Schmidt discord (HSD) [67],(24)
D
h
s
a
ρ
=
min
ρ
cq
⁡
ρ
-
ρ
c
q
2
2
,is drawing much attention due to its amenability for analytical computations when compared with most other QD measures. In the last equation, the minimization is performed over the classical-quantum states:(25)
ρ
c
q
=
∑
j
p
j
a
j
a
j
⊗
ρ
j
b
,with p
j being a probability distribution, |
a
j
〉 being an orthonormal basis for H
a, ρ
j
b being generic density operators defined in H
b, and O
2
≔
T
r
(
O
†
O
) being the Hilbert-Schmidt norm of the linear operator O, with O
† being the transpose conjugate of O.In this paper, as a basic test for the Fortran code provided to obtain coherence vectors and correlation matrices, we shall compute the following lower bound for the HSD [68]:(26)
D
h
s
a
ρ
=
∑
j
=
d
a
d
a
2
-
1
λ
j
a
,where λ
j
a are the eigenvalues, sorted in nonincreasing order, of the d
a
2
-
1
×
(
d
a
2
-
1
) matrix:(27)
Ξ
a
=
2
d
a
2
d
b
a
a
t
+
2
d
b
C
C
t
.In the above equation, t stands for the transpose of a vector or matrix. We observe that the other version of the HSD, D
h
s
b, can be obtained from the above equations simply by exchanging a and b and using C
t
C instead of C
C
t.It is interesting that, as was proved in [69], a bipartite state ρ, with polarization vectors a and b and correlation matrix C, is classical-quantum if and only if there exists a (
d
a
-
1
)-dimensional projector Π
a in the space R
d
a
2
-
1 such that(28)
Π
a
a
=
a
,
Π
a
C
=
C
.Based on this fact, an ameliorated version for the Hilbert-Schmidt quantum discord (AHSD) was proposed [69]:(29)
D
h
s
a
a
ρ
≔
min
Π
a
⁡
Υ
a
-
Π
a
Υ
a
2
2
, with the matrix Υ
a defined as(30)
Υ
a
≔
2
d
a
2
d
b
f
b
a
g
b
2
d
b
C
,where f and g are arbitrary functions of b
≡
b
2. Then, by setting f
(
b
)
=
g
(
b
)
=
P
(
ρ
b
) and using the purity,(31)
P
ρ
b
≔
T
r
ρ
b
2
∑
j
,
k
ρ
j
,
k
b
2
,to address the problem of noncontractivity of the Hilbert-Schmidt distance, the following analytical formula was presented [69]:(32)
D
h
s
a
a
ρ
=
1
P
ρ
b
∑
j
=
d
a
d
a
2
-
1
λ
j
a
=
D
h
s
a
ρ
P
ρ
b
.Thus both discord quantifiersD
h
s
a and D
h
s
a
a are, at the end of the day, obtained from the eigenvalues λ
j
a. And the computation of these eigenvalues requires the knowledge of the coherence vector a (or b) and of the correlation matrix C. These QD measures were implemented in the Fortran functions discord_hs(ssys,d
a, d
b, ρ) and discord_hsa(ssys,d
a, d
b, ρ), where ssys = ”s”, with s
=
a
,
b, specifies which version of the quantum discord is to be computed. As an example, let us use the formulas provided in this paper and the associated code to compute the HSD and AHSD of Werner states in H
a
⊗
H
b (with d
a
=
d
b): (33)
ρ
w
=
d
a
-
w
d
a
d
a
2
-
1
I
a
⊗
I
b
+
d
a
w
-
1
d
a
d
a
2
-
1
F
,where w
∈
[
-
1,1
] and F
=
∑
j
,
k
=
1
d
a
j
k
k
j
. The reduced states of ρ
w are I
s
/
d
s, whose purity is P
(
ρ
s
)
=
1
/
d
s. The results for the HSD and AHSD of ρ
w are presented in Figure 1.Figure 1
The points are the values of the ameliorated Hilbert-Schmidt quantum discord of Werner states computed numerically using (32). The lines are the corresponding values of the AHSD plotted via the analytical formula: D
h
s
a
a
(
ρ
w
)
=
d
a
D
h
s
a
(
ρ
w
)
=
(
d
a
w
-
1
)
2
/
(
(
d
a
-
1
)
(
d
a
+
1
)
2
). Due to the symmetry of ρ
w, here D
h
s
a
a
(
ρ
w
)
=
D
h
s
a
b
(
ρ
w
). In the inset the difference between the times taken by the two methods to compute the AHSD for a fixed value of d
a is shown. We see clearly here that our optimized algorithm gives an exponential speedup against the brute force calculation of Bloch’s vectors and correlation matrix.


## 4. Concluding Remarks

In this paper, we addressed the problem of computing coherence vectors and correlations matrices. We obtained formulas for these functions which make possible considerably more efficient numerical implementation when compared with the direct use of their definitions. We provided Fortran code to calculate all the quantities regarded in this paper. As a test for our formulas and code, we computed Hilbert-Schmidt quantum discords of Werner states. It is important to observe that although our focus here was in quantum information science, the tools provided can find application in other areas, such as the calculation of order parameters and correlations functions for the study of phase transitions in discrete quantum or classical systems.



---

*Source: 6892178-2016-06-29.xml*