Tomographer  v5.4
Tomographer C++ Framework Documentation
pos_semidef_util.h
Go to the documentation of this file.
1 /* This file is part of the Tomographer project, which is distributed under the
2  * terms of the MIT license.
3  *
4  * The MIT License (MIT)
5  *
6  * Copyright (c) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
7  * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, Philippe Faist
8  *
9  * Permission is hereby granted, free of charge, to any person obtaining a copy
10  * of this software and associated documentation files (the "Software"), to deal
11  * in the Software without restriction, including without limitation the rights
12  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13  * copies of the Software, and to permit persons to whom the Software is
14  * furnished to do so, subject to the following conditions:
15  *
16  * The above copyright notice and this permission notice shall be included in
17  * all copies or substantial portions of the Software.
18  *
19  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25  * SOFTWARE.
26  */
27 
28 
29 #ifndef _TOMOGRAPHER_MATHTOOLS_POS_SEMIDEF_UTIL
30 #define _TOMOGRAPHER_MATHTOOLS_POS_SEMIDEF_UTIL
31 
38 #include <Eigen/Core>
39 #include <Eigen/Eigenvalues>
40 
41 #include <tomographer/tools/cxxutil.h> // tomographer_assert()
42 
43 
44 namespace Tomographer {
45 namespace MathTools {
46 
47 
56 template<typename VectorType>
58  const typename Eigen::NumTraits<typename VectorType::Scalar>::Real tolerance = 1e-8)
59 {
60  typedef typename Eigen::NumTraits<typename VectorType::Scalar>::Real RealScalar;
61  typedef typename VectorType::Index Index;
62 
63  RealScalar trace = 0;
64  RealScalar new_excess_trace = 0;
65  Index k;
66  Index num_pos = 0;
67  for (k = 0; k < vec.size(); ++k) {
68  trace += vec(k);
69  if (vec(k) < tolerance) {
70  new_excess_trace += tolerance - vec(k);
71  } else {
72  ++num_pos;
73  }
74  }
75 
76  tomographer_assert(trace >= vec.size()*tolerance);
77 
78  RealScalar remove_from_each = new_excess_trace / num_pos;
79  for (k = 0; k < vec.size(); ++k) {
80  if (vec(k) < tolerance) {
81  vec(k) = tolerance;
82  } else {
83  vec(k) -= remove_from_each;
84  }
85  }
86 }
87 
95 template<typename MatrixType, typename Der>
96 inline MatrixType forcePosSemiDef(const Eigen::MatrixBase<Der> & rho,
97  const typename Eigen::NumTraits<typename MatrixType::Scalar>::Real tolerance = 1e-8)
98 {
99  typedef typename Eigen::SelfAdjointEigenSolver<MatrixType>::RealVectorType RealVectorType;
100 
102  rho.template cast<typename MatrixType::Scalar>()
103  );
104  auto U = eig_rho_ref.eigenvectors();
105  RealVectorType d = eig_rho_ref.eigenvalues();
106 
107  forcePosVecKeepSum<RealVectorType>(d, tolerance);
108 
109  return U * d.asDiagonal() * U.adjoint();
110 }
111 
112 
120 template<typename MatrixType>
121 inline MatrixType safeOperatorSqrt(const Eigen::Ref<const MatrixType> & A,
123  tolerance = 1e-8)
124 {
126  auto U = eig_A_ref.eigenvectors();
127  auto d = eig_A_ref.eigenvalues();
128 
129  typedef typename Eigen::SelfAdjointEigenSolver<MatrixType>::RealVectorType RealVectorType;
130 
131  forcePosVecKeepSum<RealVectorType>(d, tolerance);
132 
133  return U * d.cwiseSqrt().asDiagonal() * U.adjoint();
134 }
135 
144 template<typename MatrixType>
147  tolerance = 1e-8)
148 {
150  auto U = eig_A_ref.eigenvectors();
151  auto d = eig_A_ref.eigenvalues();
152 
153  typedef typename Eigen::SelfAdjointEigenSolver<MatrixType>::RealVectorType RealVectorType;
154 
155  forcePosVecKeepSum<RealVectorType>(d, tolerance);
156 
157  for (typename MatrixType::Index k = 0; k < d.size(); ++k) {
158  if (d(k) > tolerance) {
159  d(k) = 1 / std::sqrt(d(k));
160  }
161  }
162 
163  return U * d.asDiagonal() * U.adjoint();
164 }
165 
166 
167 
168 
169 } // MathTools
170 } // Tomographer
171 
172 
173 
174 
175 
176 #endif
Base namespace for the Tomographer project.
Definition: densellh.h:45
const RealVectorType & eigenvalues() const
MatrixType safeOperatorSqrt(const Eigen::Ref< const MatrixType > &A, const typename Eigen::NumTraits< typename MatrixType::Scalar >::Real tolerance=1e-8)
Safe version of operator square root for positive semidefinite matrices.
Some C++ utilities, with a tad of C++11 tricks.
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
void forcePosVecKeepSum(Eigen::Ref< VectorType > vec, const typename Eigen::NumTraits< typename VectorType::Scalar >::Real tolerance=1e-8)
Make sure that the given vector is numerically positive.
const MatrixType & eigenvectors() const
MatrixType safeOperatorInvSqrt(const Eigen::Ref< const MatrixType > &A, const typename Eigen::NumTraits< typename MatrixType::Scalar >::Real tolerance=1e-8)
Safe version of operator inverse square root for positive semidefinite matrices.
MatrixType forcePosSemiDef(const Eigen::MatrixBase< Der > &rho, const typename Eigen::NumTraits< typename MatrixType::Scalar >::Real tolerance=1e-8)
Make sure that rho is numerically positive semidefinite.
T sqrt(T... args)
#define tomographer_assert(...)
Assertion test macro.
Definition: cxxdefs.h:84