Quantum++  v1.0-rc2
A modern C++11 quantum computing library
entanglement.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 ENTANGLEMENT_H_
28 #define ENTANGLEMENT_H_
29 
30 namespace qpp
31 {
43 template<typename Derived>
44 dyn_col_vect<double> schmidtcoeffs(const Eigen::MatrixBase<Derived>& A,
45  const std::vector<idx>& dims)
46 {
47  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
48 
49  // EXCEPTION CHECKS
50 
51  // check zero-size
53  throw exception::ZeroSize("qpp::schmidtcoeffs()");
54  // check bi-partite
55  if (dims.size() != 2)
56  throw exception::NotBipartite("qpp::schmidtcoeffs()");
57  // check column vector
58  if (!internal::check_cvector(rA))
59  throw exception::MatrixNotCvector("qpp::schmidtcoeffs()");
60  // check matching dimensions
61  if (!internal::check_dims_match_cvect(dims, rA))
62  throw exception::DimsMismatchCvector("qpp::schmidtcoeffs()");
63  // END EXCEPTION CHECKS
64 
65  return svals(transpose(reshape(rA, dims[1], dims[0])));
66 }
67 
79 template<typename Derived>
80 dyn_col_vect<double> schmidtcoeffs(const Eigen::MatrixBase<Derived>& A,
81  idx d = 2)
82 {
83  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
84 
85  // EXCEPTION CHECKS
86 
87  // check zero size
89  throw exception::ZeroSize("qpp::schmidtcoeffs()");
90 
91  // check valid dims
92  if (d < 2)
93  throw exception::DimsInvalid("qpp::schmidtcoeffs()");
94  // END EXCEPTION CHECKS
95 
96  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
97  std::vector<idx> dims(N, d); // local dimensions vector
98 
99  return schmidtcoeffs(A, dims);
100 }
101 
110 template<typename Derived>
111 cmat schmidtA(const Eigen::MatrixBase<Derived>& A,
112  const std::vector<idx>& dims)
113 {
114  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
115 
116  // EXCEPTION CHECKS
117 
118  // check zero-size
120  throw exception::ZeroSize("qpp::schmidtA()");
121  // check bi-partite
122  if (dims.size() != 2)
123  throw exception::NotBipartite("qpp::schmidtA()");
124  // check column vector
125  if (!internal::check_cvector(rA))
126  throw exception::MatrixNotCvector("qpp::schmidtA()");
127  // check matching dimensions
128  if (!internal::check_dims_match_cvect(dims, rA))
129  throw exception::DimsMismatchCvector("qpp::schmidtA()");
130  // END EXCEPTION CHECKS
131 
132  return svdU(transpose(reshape(rA, dims[1], dims[0])));
133 }
134 
143 template<typename Derived>
144 cmat schmidtA(const Eigen::MatrixBase<Derived>& A, idx d = 2)
145 {
146  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
147 
148  // EXCEPTION CHECKS
149 
150  // check zero size
152  throw exception::ZeroSize("qpp::schmidtA()");
153 
154  // check valid dims
155  if (d < 2)
156  throw exception::DimsInvalid("qpp::schmidtA()");
157  // END EXCEPTION CHECKS
158 
159  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
160  std::vector<idx> dims(N, d); // local dimensions vector
161 
162  return schmidtA(A, dims);
163 }
164 
173 template<typename Derived>
174 cmat schmidtB(const Eigen::MatrixBase<Derived>& A,
175  const std::vector<idx>& dims)
176 {
177  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
178 
179  // EXCEPTION CHECKS
180 
181  // check zero-size
183  throw exception::ZeroSize("qpp::schmidtB()");
184  // check bi-partite
185  if (dims.size() != 2)
186  throw exception::NotBipartite("qpp::schmidtB()");
187  // check column vector
188  if (!internal::check_cvector(rA))
189  throw exception::MatrixNotCvector("qpp::schmidtB()");
190  // check matching dimensions
191  if (!internal::check_dims_match_cvect(dims, rA))
192  throw exception::DimsMismatchCvector("qpp::schmidtB()");
193  // END EXCEPTION CHECKS
194 
195  // by default returns U_B^*, we need U_B, i.e. the complex conjugate,
196  // i.e. adjoint(transpose(U_B))
197  return svdV(transpose(reshape(conjugate(rA), dims[1], dims[0])));
198 }
199 
208 template<typename Derived>
209 cmat schmidtB(const Eigen::MatrixBase<Derived>& A, idx d = 2)
210 {
211  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
212 
213  // EXCEPTION CHECKS
214 
215  // check zero size
217  throw exception::ZeroSize("qpp::schmidtB()");
218 
219  // check valid dims
220  if (d < 2)
221  throw exception::DimsInvalid("qpp::schmidtB()");
222  // END EXCEPTION CHECKS
223 
224  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
225  std::vector<idx> dims(N, d); // local dimensions vector
226 
227  return schmidtB(A, dims);
228 }
229 
242 template<typename Derived>
243 std::vector<double> schmidtprobs(const Eigen::MatrixBase<Derived>& A,
244  const std::vector<idx>& dims)
245 {
246  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
247 
248  // EXCEPTION CHECKS
249 
250  // check zero-size
252  throw exception::ZeroSize("qpp::schmidtprobs()");
253  // check bi-partite
254  if (dims.size() != 2)
255  throw exception::NotBipartite("qpp::schmidtprobs()");
256  // check column vector
257  if (!internal::check_cvector(rA))
258  throw exception::MatrixNotCvector("qpp::schmidtprobs()");
259  // check matching dimensions
260  if (!internal::check_dims_match_cvect(dims, rA))
261  throw exception::DimsMismatchCvector("qpp::schmidtprobs()");
262  // END EXCEPTION CHECKS
263 
264  std::vector<double> result;
265  dyn_col_vect<double> scf = schmidtcoeffs(rA, dims);
266  for (idx i = 0; i < static_cast<idx>(scf.rows()); ++i)
267  result.push_back(std::pow(scf(i), 2));
268 
269  return result;
270 }
271 
284 template<typename Derived>
285 std::vector<double> schmidtprobs(const Eigen::MatrixBase<Derived>& A, idx d = 2)
286 {
287  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
288 
289  // EXCEPTION CHECKS
290 
291  // check zero size
293  throw exception::ZeroSize("qpp::schmidtprobs()");
294 
295  // check valid dims
296  if (d < 2)
297  throw exception::DimsInvalid("qpp::schmidtprobs()");
298  // END EXCEPTION CHECKS
299 
300  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
301  std::vector<idx> dims(N, d); // local dimensions vector
302 
303  return schmidtprobs(A, dims);
304 }
305 
317 template<typename Derived>
318 double entanglement(const Eigen::MatrixBase<Derived>& A,
319  const std::vector<idx>& dims)
320 {
321  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
322 
323  // EXCEPTION CHECKS
324 
325  // check zero-size
327  throw exception::ZeroSize("qpp::entanglement()");
328  // check bi-partite
329  if (dims.size() != 2)
330  throw exception::NotBipartite("qpp::entanglement()");
331  // check column vector
332  if (!internal::check_cvector(rA))
333  throw exception::MatrixNotCvector("qpp::entanglement()");
334  // check matching dimensions
335  if (!internal::check_dims_match_cvect(dims, rA))
336  throw exception::DimsMismatchCvector("qpp::entanglement()");
337  // END EXCEPTION CHECKS
338 
339  return entropy(schmidtprobs(rA, dims));
340 }
341 
353 template<typename Derived>
354 double entanglement(const Eigen::MatrixBase<Derived>& A, idx d = 2)
355 {
356  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
357 
358  // EXCEPTION CHECKS
359 
360  // check zero size
362  throw exception::ZeroSize("qpp::entanglement()");
363 
364  // check valid dims
365  if (d < 2)
366  throw exception::DimsInvalid("qpp::entanglement()");
367  // END EXCEPTION CHECKS
368 
369  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
370  std::vector<idx> dims(N, d); // local dimensions vector
371 
372  return entanglement(A, dims);
373 }
374 
386 // the G-concurrence
387 template<typename Derived>
388 double gconcurrence(const Eigen::MatrixBase<Derived>& A)
389 {
390  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
391 
392  // EXCEPTION CHECKS
393 
394  // check zero-size
396  throw exception::ZeroSize("qpp::gconcurrence()");
397  // check column vector
398  if (!internal::check_cvector(rA))
399  throw exception::MatrixNotCvector("qpp::gconcurrence()");
400 
401  idx d = internal::get_dim_subsys(static_cast<idx>(rA.rows()), 2);
402 
403  // check equal local dimensions
404  if (d * d != static_cast<idx>(rA.rows()))
405  throw exception::DimsNotEqual("qpp::gconcurrence()");
406  // END EXCEPTION CHECKS
407 
408  // we compute exp(logdet()) to avoid underflow
409  return d * std::abs(std::exp(2. / d * logdet(reshape(rA, d, d))));
410 }
411 
419 template<typename Derived>
420 double negativity(const Eigen::MatrixBase<Derived>& A,
421  const std::vector<idx>& dims)
422 {
423  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
424 
425  // EXCEPTION CHECKS
426 
427  // check zero-size
429  throw exception::ZeroSize("qpp::negativity()");
430  // check bi-partite
431  if (dims.size() != 2)
432  throw exception::NotBipartite("qpp::negativity()");
433  // check square matrix vector
435  throw exception::MatrixNotSquare("qpp::negativity()");
436  // check matching dimensions
437  if (!internal::check_dims_match_mat(dims, rA))
438  throw exception::DimsMismatchMatrix("qpp::negativity()");
439  // END EXCEPTION CHECKS
440 
441  return (schatten(ptranspose(rA, {0}, dims), 1) - 1.) / 2.;
442 }
443 
451 template<typename Derived>
452 double negativity(const Eigen::MatrixBase<Derived>& A, idx d = 2)
453 {
454  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
455 
456  // EXCEPTION CHECKS
457 
458  // check zero size
460  throw exception::ZeroSize("qpp::negativity()");
461 
462  // check valid dims
463  if (d < 2)
464  throw exception::DimsInvalid("qpp::negativity()");
465  // END EXCEPTION CHECKS
466 
467  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
468  std::vector<idx> dims(N, d); // local dimensions vector
469 
470  return negativity(A, dims);
471 }
472 
480 template<typename Derived>
481 double lognegativity(const Eigen::MatrixBase<Derived>& A,
482  const std::vector<idx>& dims)
483 {
484  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
485 
486  // EXCEPTION CHECKS
487 
488  // check zero-size
490  throw exception::ZeroSize("qpp::lognegativity()");
491  // check bi-partite
492  if (dims.size() != 2)
493  throw exception::NotBipartite("qpp::lognegativity()");
494  // check square matrix vector
496  throw exception::MatrixNotSquare("qpp::lognegativity()");
497  // check matching dimensions
498  if (!internal::check_dims_match_mat(dims, rA))
499  throw exception::DimsMismatchMatrix("qpp::lognegativity()");
500  // END EXCEPTION CHECKS
501 
502  return std::log2(2 * negativity(rA, dims) + 1);
503 }
504 
512 template<typename Derived>
513 double lognegativity(const Eigen::MatrixBase<Derived>& A, idx d = 2)
514 {
515  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
516 
517  // EXCEPTION CHECKS
518 
519  // check zero size
521  throw exception::ZeroSize("qpp::lognegativity()");
522 
523  // check valid dims
524  if (d < 2)
525  throw exception::DimsInvalid("qpp::lognegativity()");
526  // END EXCEPTION CHECKS
527 
528  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
529  std::vector<idx> dims(N, d); // local dimensions vector
530 
531  return lognegativity(A, dims);
532 }
533 
540 template<typename Derived>
541 double concurrence(const Eigen::MatrixBase<Derived>& A)
542 {
543  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
544 
545  // EXCEPTION CHECKS
546 
547  // check zero-size
549  throw exception::ZeroSize("qpp::concurrence()");
550  // check square matrix vector
552  throw exception::MatrixNotSquare("qpp::concurrence()");
553  // check that the state is a 2-qubit state
554  if (rA.rows() != 4)
555  throw exception::NotQubitSubsys("qpp::concurrence()");
556  // END EXCEPTION CHECKS
557 
558  cmat sigmaY = Gates::get_instance().Y;
559  dyn_col_vect<double> lambdas =
560  evals(rA * kron(sigmaY, sigmaY) * conjugate(rA)
561  * kron(sigmaY, sigmaY)).real();
562 
563  std::vector<double> lambdas_sorted(lambdas.data(),
564  lambdas.data() + lambdas.size());
565 
566  std::sort(std::begin(lambdas_sorted), std::end(lambdas_sorted),
567  std::greater<double>());
568  std::transform(std::begin(lambdas_sorted), std::end(lambdas_sorted),
569  std::begin(lambdas_sorted), [](double elem)
570  {
571  return std::sqrt(std::abs(elem));
572  }); // chop tiny negatives
573 
574  return std::max(0.,
575  lambdas_sorted[0] - lambdas_sorted[1] - lambdas_sorted[2]
576  - lambdas_sorted[3]);
577 }
578 
579 } /* namespace qpp */
580 
581 #endif /* ENTANGLEMENT_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:129
Dimensions not equal exception.
Definition: exception.h:304
Dimension(s) mismatch matrix size exception.
Definition: exception.h:322
Derived::Scalar logdet(const Eigen::MatrixBase< Derived > &A)
Logarithm of the determinant.
Definition: functions.h:170
std::vector< double > schmidtprobs(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Schmidt probabilities of the bi-partite pure state A.
Definition: entanglement.h:243
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:122
double gconcurrence(const Eigen::MatrixBase< Derived > &A)
G-concurrence of the bi-partite pure state A.
Definition: entanglement.h:388
double lognegativity(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Logarithmic negativity of the bi-partite mixed state A.
Definition: entanglement.h:481
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
Matrix is not a column vector exception.
Definition: exception.h:168
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
Subsystems are not qubits exception.
Definition: exception.h:515
double schatten(const Eigen::MatrixBase< Derived > &A, double p)
Schatten matrix norm.
Definition: functions.h:825
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:174
Not bi-partite exception.
Definition: exception.h:532
cmat svdU(const Eigen::MatrixBase< Derived > &A)
Left singular vectors.
Definition: functions.h:494
cmat schmidtA(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Schmidt basis on Alice side.
Definition: entanglement.h:111
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:470
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:899
double entropy(const Eigen::MatrixBase< Derived > &A)
von-Neumann entropy of the density matrix A
Definition: entropies.h:39
Dimension(s) mismatch column vector size exception.
Definition: exception.h:340
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:41
dyn_col_vect< cplx > evals(const Eigen::MatrixBase< Derived > &A)
Eigenvalues.
Definition: functions.h:302
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:404
dyn_mat< typename Derived::Scalar > conjugate(const Eigen::MatrixBase< Derived > &A)
Complex conjugate.
Definition: functions.h:61
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:101
cmat schmidtB(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Schmidt basis on Bob side.
Definition: entanglement.h:174
cmat svdV(const Eigen::MatrixBase< Derived > &A)
Right singular vectors.
Definition: functions.h:519
double entanglement(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Entanglement of the bi-partite pure state A.
Definition: entanglement.h:318
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:89
static const Gates & get_instance() noexcept(std::is_nothrow_constructible< const Gates >::value)
Definition: singleton.h:89
dyn_mat< typename Derived::Scalar > ptranspose(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Partial transpose.
Definition: operations.h:1495
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:60
std::size_t idx
Non-negative integer index.
Definition: types.h:35
Matrix is not square exception.
Definition: exception.h:151
double concurrence(const Eigen::MatrixBase< Derived > &A)
Wootters concurrence of the bi-partite qubit mixed state A.
Definition: entanglement.h:541
cmat Y
Pauli Sigma-Y gate.
Definition: gates.h:45
dyn_col_vect< double > schmidtcoeffs(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Schmidt coefficients of the bi-partite pure state A.
Definition: entanglement.h:44
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1139
double negativity(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Negativity of the bi-partite mixed state A.
Definition: entanglement.h:420
Object has zero size exception.
Definition: exception.h:134