Quantum++  v1.0.0-beta2
C++11 quantum computing library
entropies.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2017 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 {
45  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
46 
47  // EXCEPTION CHECKS
48 
49  // check zero-size
51  throw Exception("qpp::entropy()", Exception::Type::ZERO_SIZE);
52 
53  // check square matrix
55  throw Exception("qpp::entropy()", Exception::Type::MATRIX_NOT_SQUARE);
56  // END EXCEPTION CHECKS
57 
58  dmat ev = svals(rA); // get the singular values
59  double result = 0;
60  for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
61  if (ev(i) != 0) // not identically zero
62  result -= ev(i) * std::log2(ev(i));
63 
64  return result;
65 }
66 
73 inline double entropy(const std::vector<double>& prob)
74 {
75  // EXCEPTION CHECKS
76 
77  // check zero-size
79  throw Exception("qpp::entropy()", Exception::Type::ZERO_SIZE);
80  // END EXCEPTION CHECKS
81 
82  double result = 0;
83  for (idx i = 0; i < prob.size(); ++i)
84  if (std::abs(prob[i]) != 0) // not identically zero
85  result -= std::abs(prob[i]) * std::log2(std::abs(prob[i]));
86 
87  return result;
88 }
89 
102 template<typename Derived>
103 double renyi(const Eigen::MatrixBase<Derived>& A, double alpha)
104 {
105  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
106 
107  // EXCEPTION CHECKS
108 
109  // check zero-size
111  throw Exception("qpp::renyi()", Exception::Type::ZERO_SIZE);
112 
113  // check square matrix
115  throw Exception("qpp::renyi()", Exception::Type::MATRIX_NOT_SQUARE);
116 
117  if (alpha < 0)
118  throw Exception("qpp::renyi()", Exception::Type::OUT_OF_RANGE);
119  // END EXCEPTION CHECKS
120 
121  if (alpha == 0) // H max
122  return std::log2(rA.rows());
123 
124  if (alpha == 1) // Shannon/von-Neumann
125  return entropy(rA);
126 
127  if (alpha == infty) // H min
128  return -std::log2(svals(rA)[0]);
129 
130  dmat sv = svals(rA); // get the singular values
131  double result = 0;
132  for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
133  if (sv(i) != 0) // not identically zero
134  result += std::pow(sv(i), alpha);
135 
136  return std::log2(result) / (1 - alpha);
137 }
138 
151 inline double renyi(const std::vector<double>& prob, double alpha)
152 {
153  // EXCEPTION CHECKS
154 
155  // check zero-size
156  if (!internal::check_nonzero_size(prob))
157  throw Exception("qpp::renyi()", Exception::Type::ZERO_SIZE);
158 
159  if (alpha < 0)
160  throw Exception("qpp::renyi()", Exception::Type::OUT_OF_RANGE);
161 
162  if (alpha == 0) // H max
163  return std::log2(prob.size());
164 
165  if (alpha == 1) // Shannon/von-Neumann
166  return entropy(prob);
167  // END EXCEPTION CHECKS
168 
169  if (alpha == infty) // H min
170  {
171  double max = 0;
172  for (idx i = 0; i < prob.size(); ++i)
173  if (std::abs(prob[i]) > max)
174  max = std::abs(prob[i]);
175 
176  return -std::log2(max);
177  }
178 
179  double result = 0;
180  for (idx i = 0; i < prob.size(); ++i)
181  if (std::abs(prob[i]) != 0) // not identically zero
182  result += std::pow(std::abs(prob[i]), alpha);
183 
184  return std::log2(result) / (1 - alpha);
185 }
186 
199 template<typename Derived>
200 double tsallis(const Eigen::MatrixBase<Derived>& A, double q)
201 {
202  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
203 
204  // EXCEPTION CHECKS
205 
206  // check zero-size
208  throw Exception("qpp::tsallis()", Exception::Type::ZERO_SIZE);
209 
210  // check square matrix
212  throw Exception("qpp::tsallis()", Exception::Type::MATRIX_NOT_SQUARE);
213 
214  if (q < 0)
215  throw Exception("qpp::tsallis()", Exception::Type::OUT_OF_RANGE);
216  // END EXCEPTION CHECKS
217 
218  if (q == 1) // Shannon/von-Neumann with base e logarithm
219  return entropy(rA) * std::log(2.);
220 
221  dmat ev = svals(rA); // get the singular values
222  double result = 0;
223  for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
224  if (ev(i) != 0) // not identically zero
225  result += std::pow(ev(i), q);
226 
227  return (result - 1) / (1 - q);
228 }
229 
242 inline double tsallis(const std::vector<double>& prob, double q)
243 {
244  // EXCEPTION CHECKS
245 
246  // check zero-size
247  if (!internal::check_nonzero_size(prob))
248  throw Exception("qpp::tsallis()", Exception::Type::ZERO_SIZE);
249 
250  if (q < 0)
251  throw Exception("qpp::tsallis()", Exception::Type::OUT_OF_RANGE);
252  // END EXCEPTION CHECKS
253 
254  if (q == 1) // Shannon/von-Neumann with base e logarithm
255  return entropy(prob) * std::log(2.);
256 
257  double result = 0;
258  for (idx i = 0; i < prob.size(); ++i)
259  if (std::abs(prob[i]) != 0) // not identically zero
260  result += std::pow(std::abs(prob[i]), q);
261 
262  return (result - 1) / (1 - q);
263 }
264 
274 template<typename Derived>
275 double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
276  const std::vector<idx>& subsysA,
277  const std::vector<idx>& subsysB,
278  const std::vector<idx>& dims)
279 {
280  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
281 
282  // EXCEPTION CHECKS
283 
284  // check zero-size
286  throw Exception("qpp::qmutualinfo()", Exception::Type::ZERO_SIZE);
287 
288  // check that dims is a valid dimension vector
289  if (!internal::check_dims(dims))
290  throw Exception("qpp::qmutualinfo()", Exception::Type::DIMS_INVALID);
291 
292  // check square matrix
294  throw Exception("qpp::qmutualinfo()",
296 
297  // check that dims match the dimension of A
298  if (!internal::check_dims_match_mat(dims, rA))
299  throw Exception("qpp::qmutualinfo()",
301 
302  // check that subsys are valid
303  if (!internal::check_subsys_match_dims(subsysA, dims)
304  || !internal::check_subsys_match_dims(subsysB, dims))
305  throw Exception("qpp::qmutualinfo()",
307  // END EXCEPTION CHECKS
308 
309  // The full system indexes {0,1,...,n-1}
310  std::vector<idx> full_system(dims.size());
311  std::iota(std::begin(full_system), std::end(full_system), 0);
312 
313  // Sorted input subsystems
314  std::vector<idx> subsysAsorted{subsysA};
315  std::vector<idx> subsysBsorted{subsysB};
316 
317  // sort the input subsystems (as needed by std::set_union)
318  std::sort(std::begin(subsysAsorted), std::end(subsysAsorted));
319  std::sort(std::begin(subsysBsorted), std::end(subsysBsorted));
320 
321  // construct the union of A and B
322  std::vector<idx> subsysAB;
323  std::set_union(std::begin(subsysAsorted), std::end(subsysAsorted),
324  std::begin(subsysBsorted), std::end(subsysBsorted),
325  std::back_inserter(subsysAB));
326 
327  // construct the complements
328  std::vector<idx> subsysA_bar = complement(subsysA, dims.size());
329  std::vector<idx> subsysB_bar = complement(subsysB, dims.size());;
330  std::vector<idx> subsysAB_bar = complement(subsysAB, dims.size());
331 
332  cmat rhoA = ptrace(rA, subsysA_bar, dims);
333  cmat rhoB = ptrace(rA, subsysB_bar, dims);
334  cmat rhoAB = ptrace(rA, subsysAB_bar, dims);
335 
336  return entropy(rhoA) + entropy(rhoB) - entropy(rhoAB);
337 }
338 
348 template<typename Derived>
349 double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
350  const std::vector<idx>& subsysA,
351  const std::vector<idx>& subsysB,
352  idx d = 2)
353 {
354  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
355 
356  // EXCEPTION CHECKS
357 
358  // check zero-size
360  throw Exception("qpp::qmutualinfo()", Exception::Type::ZERO_SIZE);
361 
362  // check valid dims
363  if (d == 0)
364  throw Exception("qpp::qmutualinfo()", Exception::Type::DIMS_INVALID);
365  // END EXCEPTION CHECKS
366 
367  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
368  std::vector<idx> dims(N, d); // local dimensions vector
369 
370  return qmutualinfo(rA, subsysA, subsysB, dims);
371 }
372 
373 } /* namespace qpp */
374 
375 #endif /* ENTROPY_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:127
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:275
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:198
double tsallis(const Eigen::MatrixBase< Derived > &A, double q)
Tsallis- entropy of the density matrix A, for .
Definition: entropies.h:200
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:66
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:78
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:156
Quantum++ main namespace.
Definition: codes.h:30
idx get_num_subsys(idx sz, idx d)
Definition: util.h:370
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1814
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:470
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:103
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:1280
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:88
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:61
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:140
std::size_t idx
Non-negative integer index.
Definition: types.h:36