Quantum++  v1.0-rc4
A modern C++11 quantum computing library
entropies.h
Go to the documentation of this file.
1 /*
2  * This file is part of Quantum++.
3  *
4  * MIT License
5  *
6  * Copyright (c) 2013 - 2018 Vlad Gheorghiu (vgheorgh@gmail.com)
7  *
8  * Permission is hereby granted, free of charge, to any person obtaining a copy
9  * of this software and associated documentation files (the "Software"), to deal
10  * in the Software without restriction, including without limitation the rights
11  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12  * copies of the Software, and to permit persons to whom the Software is
13  * furnished to do so, subject to the following conditions:
14  *
15  * The above copyright notice and this permission notice shall be included in
16  * all copies or substantial portions of the Software.
17  *
18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24  * SOFTWARE.
25  */
26 
32 #ifndef ENTROPY_H_
33 #define ENTROPY_H_
34 
35 namespace qpp {
42 template <typename Derived>
43 double entropy(const Eigen::MatrixBase<Derived>& A) {
44  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
45 
46  // EXCEPTION CHECKS
47 
48  // check zero-size
50  throw exception::ZeroSize("qpp::entropy()");
51 
52  // check square matrix
54  throw exception::MatrixNotSquare("qpp::entropy()");
55  // END EXCEPTION CHECKS
56 
57  dmat ev = svals(rA); // get the singular values
58  double result = 0;
59  for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
60  if (ev(i) != 0) // not identically zero
61  result -= ev(i) * std::log2(ev(i));
62 
63  return result;
64 }
65 
72 inline double entropy(const std::vector<double>& prob) {
73  // EXCEPTION CHECKS
74 
75  // check zero-size
77  throw exception::ZeroSize("qpp::entropy()");
78  // END EXCEPTION CHECKS
79 
80  double result = 0;
81  for (idx i = 0; i < prob.size(); ++i)
82  if (std::abs(prob[i]) != 0) // not identically zero
83  result -= std::abs(prob[i]) * std::log2(std::abs(prob[i]));
84 
85  return result;
86 }
87 
100 template <typename Derived>
101 double renyi(const Eigen::MatrixBase<Derived>& A, double alpha) {
102  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
103 
104  // EXCEPTION CHECKS
105 
106  // check zero-size
108  throw exception::ZeroSize("qpp::renyi()");
109 
110  // check square matrix
112  throw exception::MatrixNotSquare("qpp::renyi()");
113 
114  if (alpha < 0)
115  throw exception::OutOfRange("qpp::renyi()");
116  // END EXCEPTION CHECKS
117 
118  if (alpha == 0) // H max
119  return std::log2(rA.rows());
120 
121  if (alpha == 1) // Shannon/von-Neumann
122  return entropy(rA);
123 
124  if (alpha == infty) // H min
125  return -std::log2(svals(rA)[0]);
126 
127  dmat sv = svals(rA); // get the singular values
128  double result = 0;
129  for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
130  if (sv(i) != 0) // not identically zero
131  result += std::pow(sv(i), alpha);
132 
133  return std::log2(result) / (1 - alpha);
134 }
135 
148 inline double renyi(const std::vector<double>& prob, double alpha) {
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  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
198 
199  // EXCEPTION CHECKS
200 
201  // check zero-size
203  throw exception::ZeroSize("qpp::tsallis()");
204 
205  // check square matrix
207  throw exception::MatrixNotSquare("qpp::tsallis()");
208 
209  if (q < 0)
210  throw exception::OutOfRange("qpp::tsallis()");
211  // END EXCEPTION CHECKS
212 
213  if (q == 1) // Shannon/von-Neumann with base e logarithm
214  return entropy(rA) * std::log(2);
215 
216  dmat ev = svals(rA); // get the singular values
217  double result = 0;
218  for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
219  if (ev(i) != 0) // not identically zero
220  result += std::pow(ev(i), q);
221 
222  return (result - 1) / (1 - q);
223 }
224 
237 inline double tsallis(const std::vector<double>& prob, double q) {
238  // EXCEPTION CHECKS
239 
240  // check zero-size
241  if (!internal::check_nonzero_size(prob))
242  throw exception::ZeroSize("qpp::tsallis()");
243 
244  if (q < 0)
245  throw exception::OutOfRange("qpp::tsallis()");
246  // END EXCEPTION CHECKS
247 
248  if (q == 1) // Shannon/von-Neumann with base e logarithm
249  return entropy(prob) * std::log(2);
250 
251  double result = 0;
252  for (idx i = 0; i < prob.size(); ++i)
253  if (std::abs(prob[i]) != 0) // not identically zero
254  result += std::pow(std::abs(prob[i]), q);
255 
256  return (result - 1) / (1 - q);
257 }
258 
268 template <typename Derived>
269 double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
270  const std::vector<idx>& subsysA,
271  const std::vector<idx>& subsysB,
272  const std::vector<idx>& dims) {
273  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
274 
275  // EXCEPTION CHECKS
276 
277  // check zero-size
279  throw exception::ZeroSize("qpp::qmutualinfo()");
280 
281  // check that dims is a valid dimension vector
282  if (!internal::check_dims(dims))
283  throw exception::DimsInvalid("qpp::qmutualinfo()");
284 
285  // check square matrix
287  throw exception::MatrixNotSquare("qpp::qmutualinfo()");
288 
289  // check that dims match the dimension of A
290  if (!internal::check_dims_match_mat(dims, rA))
291  throw exception::DimsMismatchMatrix("qpp::qmutualinfo()");
292 
293  // check that subsys are valid
294  if (!internal::check_subsys_match_dims(subsysA, dims) ||
295  !internal::check_subsys_match_dims(subsysB, dims))
296  throw exception::SubsysMismatchDims("qpp::qmutualinfo()");
297  // END EXCEPTION CHECKS
298 
299  // The full system indexes {0,1,...,n-1}
300  std::vector<idx> full_system(dims.size());
301  std::iota(std::begin(full_system), std::end(full_system), 0);
302 
303  // Sorted input subsystems
304  std::vector<idx> subsysAsorted{subsysA};
305  std::vector<idx> subsysBsorted{subsysB};
306 
307  // sort the input subsystems (as needed by std::set_union)
308  std::sort(std::begin(subsysAsorted), std::end(subsysAsorted));
309  std::sort(std::begin(subsysBsorted), std::end(subsysBsorted));
310 
311  // construct the union of A and B
312  std::vector<idx> subsysAB;
313  std::set_union(std::begin(subsysAsorted), std::end(subsysAsorted),
314  std::begin(subsysBsorted), std::end(subsysBsorted),
315  std::back_inserter(subsysAB));
316 
317  // construct the complements
318  std::vector<idx> subsysA_bar = complement(subsysA, dims.size());
319  std::vector<idx> subsysB_bar = complement(subsysB, dims.size());
320  ;
321  std::vector<idx> subsysAB_bar = complement(subsysAB, dims.size());
322 
323  cmat rhoA = ptrace(rA, subsysA_bar, dims);
324  cmat rhoB = ptrace(rA, subsysB_bar, dims);
325  cmat rhoAB = ptrace(rA, subsysAB_bar, dims);
326 
327  return entropy(rhoA) + entropy(rhoB) - entropy(rhoAB);
328 }
329 
339 template <typename Derived>
340 double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
341  const std::vector<idx>& subsysA,
342  const std::vector<idx>& subsysB, idx d = 2) {
343  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
344 
345  // EXCEPTION CHECKS
346 
347  // check zero-size
349  throw exception::ZeroSize("qpp::qmutualinfo()");
350 
351  // check valid dims
352  if (d < 2)
353  throw exception::DimsInvalid("qpp::qmutualinfo()");
354  // END EXCEPTION CHECKS
355 
356  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
357  std::vector<idx> dims(N, d); // local dimensions vector
358 
359  return qmutualinfo(rA, subsysA, subsysB, dims);
360 }
361 
362 } /* namespace qpp */
363 
364 #endif /* ENTROPY_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
Dimension(s) mismatch matrix size exception.
Definition: exception.h:298
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:269
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:210
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:69
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:81
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:150
Subsystems mismatch dimensions exception.
Definition: exception.h:363
Quantum++ main namespace.
Definition: codes.h:35
idx get_num_subsys(idx sz, idx d)
Definition: util.h:366
Invalid dimension(s) exception.
Definition: exception.h:267
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1733
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:455
double entropy(const Eigen::MatrixBase< Derived > &A)
von-Neumann entropy of the density matrix A
Definition: entropies.h:43
double renyi(const Eigen::MatrixBase< Derived > &A, double alpha)
Renyi- entropy of the density matrix A, for .
Definition: entropies.h:101
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:1176
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:87
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
Parameter out of range exception.
Definition: exception.h:513
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
std::size_t idx
Non-negative integer index.
Definition: types.h:39
Matrix is not square exception.
Definition: exception.h:147
Object has zero size exception.
Definition: exception.h:132