Quantum++  v1.0-rc2
A modern 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 namespace qpp
31 {
38 template<typename Derived>
39 double entropy(const Eigen::MatrixBase<Derived>& A)
40 {
41  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
42 
43  // EXCEPTION CHECKS
44 
45  // check zero-size
47  throw exception::ZeroSize("qpp::entropy()");
48 
49  // check square matrix
51  throw exception::MatrixNotSquare("qpp::entropy()");
52  // END EXCEPTION CHECKS
53 
54  dmat ev = svals(rA); // get the singular values
55  double result = 0;
56  for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
57  if (ev(i) != 0) // not identically zero
58  result -= ev(i) * std::log2(ev(i));
59 
60  return result;
61 }
62 
69 inline double entropy(const std::vector<double>& prob)
70 {
71  // EXCEPTION CHECKS
72 
73  // check zero-size
75  throw exception::ZeroSize("qpp::entropy()");
76  // END EXCEPTION CHECKS
77 
78  double result = 0;
79  for (idx i = 0; i < prob.size(); ++i)
80  if (std::abs(prob[i]) != 0) // not identically zero
81  result -= std::abs(prob[i]) * std::log2(std::abs(prob[i]));
82 
83  return result;
84 }
85 
98 template<typename Derived>
99 double renyi(const Eigen::MatrixBase<Derived>& A, double alpha)
100 {
101  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
102 
103  // EXCEPTION CHECKS
104 
105  // check zero-size
107  throw exception::ZeroSize("qpp::renyi()");
108 
109  // check square matrix
111  throw exception::MatrixNotSquare("qpp::renyi()");
112 
113  if (alpha < 0)
114  throw exception::OutOfRange("qpp::renyi()");
115  // END EXCEPTION CHECKS
116 
117  if (alpha == 0) // H max
118  return std::log2(rA.rows());
119 
120  if (alpha == 1) // Shannon/von-Neumann
121  return entropy(rA);
122 
123  if (alpha == infty) // H min
124  return -std::log2(svals(rA)[0]);
125 
126  dmat sv = svals(rA); // get the singular values
127  double result = 0;
128  for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
129  if (sv(i) != 0) // not identically zero
130  result += std::pow(sv(i), alpha);
131 
132  return std::log2(result) / (1 - alpha);
133 }
134 
147 inline double renyi(const std::vector<double>& prob, double alpha)
148 {
149  // EXCEPTION CHECKS
150 
151  // check zero-size
152  if (!internal::check_nonzero_size(prob))
153  throw exception::ZeroSize("qpp::renyi()");
154 
155  if (alpha < 0)
156  throw exception::OutOfRange("qpp::renyi()");
157 
158  if (alpha == 0) // H max
159  return std::log2(prob.size());
160 
161  if (alpha == 1) // Shannon/von-Neumann
162  return entropy(prob);
163  // END EXCEPTION CHECKS
164 
165  if (alpha == infty) // H min
166  {
167  double max = 0;
168  for (idx i = 0; i < prob.size(); ++i)
169  if (std::abs(prob[i]) > max)
170  max = std::abs(prob[i]);
171 
172  return -std::log2(max);
173  }
174 
175  double result = 0;
176  for (idx i = 0; i < prob.size(); ++i)
177  if (std::abs(prob[i]) != 0) // not identically zero
178  result += std::pow(std::abs(prob[i]), alpha);
179 
180  return std::log2(result) / (1 - alpha);
181 }
182 
195 template<typename Derived>
196 double tsallis(const Eigen::MatrixBase<Derived>& A, double q)
197 {
198  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
199 
200  // EXCEPTION CHECKS
201 
202  // check zero-size
204  throw exception::ZeroSize("qpp::tsallis()");
205 
206  // check square matrix
208  throw exception::MatrixNotSquare("qpp::tsallis()");
209 
210  if (q < 0)
211  throw exception::OutOfRange("qpp::tsallis()");
212  // END EXCEPTION CHECKS
213 
214  if (q == 1) // Shannon/von-Neumann with base e logarithm
215  return entropy(rA) * std::log(2);
216 
217  dmat ev = svals(rA); // get the singular values
218  double result = 0;
219  for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
220  if (ev(i) != 0) // not identically zero
221  result += std::pow(ev(i), q);
222 
223  return (result - 1) / (1 - q);
224 }
225 
238 inline double tsallis(const std::vector<double>& prob, double q)
239 {
240  // EXCEPTION CHECKS
241 
242  // check zero-size
243  if (!internal::check_nonzero_size(prob))
244  throw exception::ZeroSize("qpp::tsallis()");
245 
246  if (q < 0)
247  throw exception::OutOfRange("qpp::tsallis()");
248  // END EXCEPTION CHECKS
249 
250  if (q == 1) // Shannon/von-Neumann with base e logarithm
251  return entropy(prob) * std::log(2);
252 
253  double result = 0;
254  for (idx i = 0; i < prob.size(); ++i)
255  if (std::abs(prob[i]) != 0) // not identically zero
256  result += std::pow(std::abs(prob[i]), q);
257 
258  return (result - 1) / (1 - q);
259 }
260 
270 template<typename Derived>
271 double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
272  const std::vector<idx>& subsysA,
273  const std::vector<idx>& subsysB,
274  const std::vector<idx>& dims)
275 {
276  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
277 
278  // EXCEPTION CHECKS
279 
280  // check zero-size
282  throw exception::ZeroSize("qpp::qmutualinfo()");
283 
284  // check that dims is a valid dimension vector
285  if (!internal::check_dims(dims))
286  throw exception::DimsInvalid("qpp::qmutualinfo()");
287 
288  // check square matrix
290  throw exception::MatrixNotSquare("qpp::qmutualinfo()");
291 
292  // check that dims match the dimension of A
293  if (!internal::check_dims_match_mat(dims, rA))
294  throw exception::DimsMismatchMatrix("qpp::qmutualinfo()");
295 
296  // check that subsys are valid
297  if (!internal::check_subsys_match_dims(subsysA, dims)
298  || !internal::check_subsys_match_dims(subsysB, dims))
299  throw exception::SubsysMismatchDims("qpp::qmutualinfo()");
300  // END EXCEPTION CHECKS
301 
302  // The full system indexes {0,1,...,n-1}
303  std::vector<idx> full_system(dims.size());
304  std::iota(std::begin(full_system), std::end(full_system), 0);
305 
306  // Sorted input subsystems
307  std::vector<idx> subsysAsorted{subsysA};
308  std::vector<idx> subsysBsorted{subsysB};
309 
310  // sort the input subsystems (as needed by std::set_union)
311  std::sort(std::begin(subsysAsorted), std::end(subsysAsorted));
312  std::sort(std::begin(subsysBsorted), std::end(subsysBsorted));
313 
314  // construct the union of A and B
315  std::vector<idx> subsysAB;
316  std::set_union(std::begin(subsysAsorted), std::end(subsysAsorted),
317  std::begin(subsysBsorted), std::end(subsysBsorted),
318  std::back_inserter(subsysAB));
319 
320  // construct the complements
321  std::vector<idx> subsysA_bar = complement(subsysA, dims.size());
322  std::vector<idx> subsysB_bar = complement(subsysB, dims.size());;
323  std::vector<idx> subsysAB_bar = complement(subsysAB, dims.size());
324 
325  cmat rhoA = ptrace(rA, subsysA_bar, dims);
326  cmat rhoB = ptrace(rA, subsysB_bar, dims);
327  cmat rhoAB = ptrace(rA, subsysAB_bar, dims);
328 
329  return entropy(rhoA) + entropy(rhoB) - entropy(rhoAB);
330 }
331 
341 template<typename Derived>
342 double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
343  const std::vector<idx>& subsysA,
344  const std::vector<idx>& subsysB,
345  idx d = 2)
346 {
347  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
348 
349  // EXCEPTION CHECKS
350 
351  // check zero-size
353  throw exception::ZeroSize("qpp::qmutualinfo()");
354 
355  // check valid dims
356  if (d < 2)
357  throw exception::DimsInvalid("qpp::qmutualinfo()");
358  // END EXCEPTION CHECKS
359 
360  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
361  std::vector<idx> dims(N, d); // local dimensions vector
362 
363  return qmutualinfo(rA, subsysA, subsysB, dims);
364 }
365 
366 } /* namespace qpp */
367 
368 #endif /* ENTROPY_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:129
Dimension(s) mismatch matrix size exception.
Definition: exception.h:322
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:271
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:221
double tsallis(const Eigen::MatrixBase< Derived > &A, double q)
Tsallis- entropy of the density matrix A, for .
Definition: entropies.h:196
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:65
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:77
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:158
Subsystems mismatch dimensions exception.
Definition: exception.h:395
Quantum++ main namespace.
Definition: codes.h:30
idx get_num_subsys(idx sz, idx d)
Definition: util.h:391
Invalid dimension(s) exception.
Definition: exception.h:287
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1811
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:39
double renyi(const Eigen::MatrixBase< Derived > &A, double alpha)
Renyi- entropy of the density matrix A, for .
Definition: entropies.h:99
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:1241
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:85
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:101
Parameter out of range exception.
Definition: exception.h:567
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:60
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:142
std::size_t idx
Non-negative integer index.
Definition: types.h:35
Matrix is not square exception.
Definition: exception.h:151
Object has zero size exception.
Definition: exception.h:134