Quantum++  v0.8
C++11 quantum computing library
entropies.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2015 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef ENTROPY_H_
28 #define ENTROPY_H_
29 
30 // various entropies, assume as input either
31 // a normalized Hermitian matrix or a probability vector
32 
33 namespace qpp
34 {
35 
42 template<typename Derived>
43 double entropy(const Eigen::MatrixBase<Derived>& A)
44 {
46 
47  // check zero-size
49  throw Exception("qpp::entropy()", Exception::Type::ZERO_SIZE);
50 
51  // check square matrix
53  throw Exception("qpp::entropy()", Exception::Type::MATRIX_NOT_SQUARE);
54 
55  dmat ev = svals(rA); // get the singular values
56  double result = 0;
57  for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
58  if (ev(i) != 0) // not identically zero
59  result -= ev(i) * std::log2(ev(i));
60 
61  return result;
62 }
63 
70 inline double entropy(const std::vector <double>& prob)
71 {
72  // check zero-size
74  throw Exception("qpp::entropy()", Exception::Type::ZERO_SIZE);
75 
76  double result = 0;
77  for (idx i = 0; i < prob.size(); ++i)
78  if (std::abs(prob[i]) != 0) // not identically zero
79  result -= std::abs(prob[i]) * std::log2(std::abs(prob[i]));
80 
81  return result;
82 }
83 
96 template<typename Derived>
97 double renyi(const Eigen::MatrixBase<Derived>& A, double alpha)
98 {
100 
101  // check zero-size
103  throw Exception("qpp::renyi()", Exception::Type::ZERO_SIZE);
104 
105  // check square matrix
107  throw Exception("qpp::renyi()", Exception::Type::MATRIX_NOT_SQUARE);
108 
109  if (alpha < 0)
110  throw Exception("qpp::renyi()", Exception::Type::OUT_OF_RANGE);
111 
112  if (alpha == 0) // H max
113  return std::log2(static_cast<double>( rA.rows()));
114 
115  if (alpha == 1) // Shannon/von-Neumann
116  return entropy(rA);
117 
118  if (alpha == infty) // H min
119  return -std::log2(svals(rA)[0]);
120 
121  dmat sv = svals(rA); // get the singular values
122  double result = 0;
123  for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
124  if (sv(i) != 0) // not identically zero
125  result += std::pow(sv(i), alpha);
126 
127  return std::log2(result) / (1 - alpha);
128 }
129 
142 inline double renyi(const std::vector <double>& prob, double alpha)
143 {
144  // check zero-size
146  throw Exception("qpp::renyi()", Exception::Type::ZERO_SIZE);
147 
148  if (alpha < 0)
149  throw Exception("qpp::renyi()", Exception::Type::OUT_OF_RANGE);
150 
151  if (alpha == 0) // H max
152  return std::log2(static_cast<double>( prob.size()));
153 
154  if (alpha == 1) // Shannon/von-Neumann
155  return entropy(prob);
156 
157  if (alpha == infty) // H min
158  {
159  double max = 0;
160  for (idx i = 0; i < prob.size(); ++i)
161  if (std::abs(prob[i]) > max)
162  max = std::abs(prob[i]);
163 
164  return -std::log2(max);
165  }
166 
167  double result = 0;
168  for (idx i = 0; i < prob.size(); ++i)
169  if (std::abs(prob[i]) != 0) // not identically zero
170  result += std::pow(std::abs(prob[i]), alpha);
171 
172  return std::log2(result) / (1 - alpha);
173 }
174 
187 template<typename Derived>
188 double tsallis(const Eigen::MatrixBase<Derived>& A, double q)
189 {
190  const dyn_mat<typename Derived::Scalar>& rA = A;
191 
192  // check zero-size
194  throw Exception("qpp::tsallis()", Exception::Type::ZERO_SIZE);
195 
196  // check square matrix
198  throw Exception("qpp::tsallis()", Exception::Type::MATRIX_NOT_SQUARE);
199 
200  if (q < 0)
201  throw Exception("qpp::tsallis()", Exception::Type::OUT_OF_RANGE);
202 
203  if (q == 1) // Shannon/von-Neumann with base e logarithm
204  return entropy(rA) * std::log(2.);
205 
206  dmat ev = svals(rA); // get the singular values
207  double result = 0;
208  for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
209  if (ev(i) != 0) // not identically zero
210  result += std::pow(ev(i), q);
211 
212  return (result - 1) / (1 - q);
213 }
214 
227 inline double tsallis(const std::vector <double>& prob, double q)
228 {
229  // check zero-size
231  throw Exception("qpp::tsallis()", Exception::Type::ZERO_SIZE);
232 
233  if (q < 0)
234  throw Exception("qpp::tsallis()", Exception::Type::OUT_OF_RANGE);
235 
236  if (q == 1) // Shannon/von-Neumann with base e logarithm
237  return entropy(prob) * std::log(2.);
238 
239  double result = 0;
240  for (idx i = 0; i < prob.size(); ++i)
241  if (std::abs(prob[i]) != 0) // not identically zero
242  result += std::pow(std::abs(prob[i]), q);
243 
244  return (result - 1) / (1 - q);
245 }
246 
256 template<typename Derived>
257 double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
258  const std::vector <idx>& subsysA,
259  const std::vector <idx>& subsysB,
260  const std::vector <idx>& dims)
261 {
262  const dyn_mat<typename Derived::Scalar>& rA = A;
263 
264  // error checks
265 
266  // check zero-size
268  throw Exception("qpp::qmutualinfo()", Exception::Type::ZERO_SIZE);
269 
270  // check that dims is a valid dimension vector
271  if (!internal::_check_dims(dims))
272  throw Exception("qpp::qmutualinfo()", Exception::Type::DIMS_INVALID);
273 
274  // check square matrix
276  throw Exception("qpp::qmutualinfo()",
278 
279  // check that dims match the dimension of A
280  if (!internal::_check_dims_match_mat(dims, rA))
281  throw Exception("qpp::qmutualinfo()",
283 
284  // check that subsys are valid
285  if (!internal::_check_subsys_match_dims(subsysA, dims)
286  || !internal::_check_subsys_match_dims(subsysB, dims))
287  throw Exception("qpp::qmutualinfo()",
289 
290  // The full system indexes {0,1,...,n-1}
291  std::vector <idx> full_system(dims.size());
292  std::iota(std::begin(full_system), std::end(full_system), 0);
293 
294  // Sorted input subsystems
295  std::vector <idx> subsysAsorted{subsysA};
296  std::vector <idx> subsysBsorted{subsysB};
297 
298  // sort the input subsystems (as needed by std::set_union)
299  std::sort(std::begin(subsysAsorted), std::end(subsysAsorted));
300  std::sort(std::begin(subsysBsorted), std::end(subsysBsorted));
301 
302  // construct the union of A and B
303  std::vector <idx> subsysAB;
304  std::set_union(std::begin(subsysAsorted), std::end(subsysAsorted),
305  std::begin(subsysBsorted), std::end(subsysBsorted),
306  std::back_inserter(subsysAB));
307 
308  // construct the complements
309  std::vector <idx> subsysAbar = complement(subsysA, dims.size());
310  std::vector <idx> subsysBbar = complement(subsysB, dims.size());;
311  std::vector <idx> subsysABbar = complement(subsysAB, dims.size());
312 
313  cmat rhoA = ptrace(rA, subsysAbar, dims);
314  cmat rhoB = ptrace(rA, subsysBbar, dims);
315  cmat rhoAB = ptrace(rA, subsysABbar, dims);
316 
317  return entropy(rhoA) + entropy(rhoB) - entropy(rhoAB);
318 }
319 
329 template<typename Derived>
330 double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
331  const std::vector <idx>& subsysA,
332  const std::vector <idx>& subsysB,
333  idx d = 2)
334 {
335  const dyn_mat<typename Derived::Scalar>& rA = A;
336 
337  // error checks
338 
339  // check zero-size
341  throw Exception("qpp::qmutualinfo()", Exception::Type::ZERO_SIZE);
342 
343  idx n =
344  static_cast<idx>(std::llround(std::log2(rA.rows()) /
345  std::log2(d)));
346  std::vector <idx> dims(n, d); // local dimensions vector
347 
348  return qmutualinfo(rA, subsysA, subsysB, dims);
349 }
350 
351 } /* namespace qpp */
352 
353 #endif /* ENTROPY_H_ */
double qmutualinfo(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &subsysA, const std::vector< idx > &subsysB, const std::vector< idx > &dims)
Quantum mutual information between 2 subsystems of a composite system.
Definition: entropies.h:257
double tsallis(const Eigen::MatrixBase< Derived > &A, double q)
Tsallis- entropy of the density matrix A, for .
Definition: entropies.h:188
bool _check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:183
bool _check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:135
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:72
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:84
Quantum++ main namespace.
Definition: codes.h:30
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:77
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1653
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:426
double entropy(const Eigen::MatrixBase< Derived > &A)
von-Neumann entropy of the density matrix A
Definition: entropies.h:43
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
double renyi(const Eigen::MatrixBase< Derived > &A, double alpha)
Renyi- entropy of the density matrix A, for .
Definition: entropies.h:97
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:113
dyn_mat< typename Derived::Scalar > ptrace(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1079
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:96
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:67
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:119