Quantum++  v0.8.8
C++11 quantum computing library
entropies.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2016 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
114  if ( !internal::_check_square_mat(rA))
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(static_cast<double>( 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(static_cast<double>( 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
211  if ( !internal::_check_square_mat(rA))
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
293  if ( !internal::_check_square_mat(rA))
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> subsysAbar = complement(subsysA, dims.size());
329  std::vector<idx> subsysBbar = complement(subsysB, dims.size());;
330  std::vector<idx> subsysABbar = complement(subsysAB, dims.size());
331 
332  cmat rhoA = ptrace(rA, subsysAbar, dims);
333  cmat rhoB = ptrace(rA, subsysBbar, dims);
334  cmat rhoAB = ptrace(rA, subsysABbar, 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 =
368  static_cast<idx>(std::llround(std::log2(rA.rows()) /
369  std::log2(d)));
370  std::vector<idx> dims(n, d); // local dimensions vector
371 
372  return qmutualinfo(rA, subsysA, subsysB, dims);
373 }
374 
375 } /* namespace qpp */
376 
377 #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:275
double tsallis(const Eigen::MatrixBase< Derived > &A, double q)
Tsallis- entropy of the density matrix A, for .
Definition: entropies.h:200
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:141
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:71
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:83
Quantum++ main namespace.
Definition: codes.h:30
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:84
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1818
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:475
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
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:112
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:1218
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:88
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:66
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:125