Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
decomposition.hpp
1 /*
2  This file is part of the Kernel Quantum Probability library (KQP).
3 
4  KQP is free software: you can redistribute it and/or modify
5  it under the terms of the GNU General Public License as published by
6  the Free Software Foundation, either version 3 of the License, or
7  (at your option) any later version.
8 
9  KQP is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with KQP. If not, see <http://www.gnu.org/licenses/>.
16  */
17 
18 #ifndef _KQP_DECOMPOSITION_H
19 #define _KQP_DECOMPOSITION_H
20 
21 #include <utility>
22 #include <kqp/feature_matrix.hpp>
23 #include <kqp/evd_utils.hpp>
24 
25 
26 namespace kqp
27 {
28 
34 template<typename Scalar>
36 {
37 
38 public:
39  KQP_SCALAR_TYPEDEFS(Scalar);
40 
42  FSpaceCPtr fs;
43 
45  FMatrixPtr mX;
46 
48  ScalarAltMatrix mY;
49 
51  RealAltVector mD;
52 
55 
58 
60  Index updateCount;
61 
63  Decomposition() : orthonormal(true), m_squareRoot(false) {}
64 
66  Decomposition(const FSpaceCPtr &fs) : fs(fs), mX(fs->newMatrix()), orthonormal(true), m_squareRoot(false) {}
67 
69  Decomposition(const FSpaceCPtr &fs, const FMatrixCPtr &mX, const ScalarAltMatrix &mY, const RealAltVector &mD, bool orthonormal, bool squareRoot)
70  : fs(fs), mX(mX->copy()), mY(mY), mD(mD), orthonormal(orthonormal), m_squareRoot(squareRoot), updateCount(0) {}
71 
72 #ifndef SWIG
73 
74  Decomposition(FSpaceCPtr && fs, FMatrixPtr && mX, const ScalarAltMatrix && mY, const RealAltVector && mD, bool orthonormal, bool squareRoot)
75  : fs(fs), mX(mX), mY(mY), mD(mD), orthonormal(orthonormal), m_squareRoot(squareRoot), updateCount(0)
76  {
77  }
78 
81  {
82  *this = std::move(other);
83  }
84 
87  {
88  take(other);
89  return *this;
90  }
91 #endif
92  void take(Decomposition &other)
93  {
94  fs = other.fs;
95  mX = other.mX;
96  mY.swap(other.mY);
97  mD.swap(other.mD);
98  std::swap(orthonormal, other.orthonormal);
99  std::swap(updateCount, other.updateCount);
100  std::swap(m_squareRoot, other.m_squareRoot);
101  }
102 
103 
106  {
107  *this = other;
108  }
109 
112  {
113  fs = other.fs;
114  mX = other.mX->copy();
115  mY = other.mY;
116  mD = other.mD;
117  orthonormal = other.orthonormal;
118  updateCount = other.updateCount;
119  m_squareRoot = other.m_squareRoot;
120  return *this;
121  }
122 
123 
127  ScalarMatrix k(const Decomposition &other) const
128  {
129  return fs->k(mX, mY, mD, other.mX, other.mY, other.mD);
130  }
131 
133  bool check() const {
134  return mX->size() == mY.rows() && mY.cols() == mD.rows();
135  }
136 
137  Index rank() const {
138  return mD.size();
139  }
140 
141  Index preimagesCount() const {
142  return mX->size();
143  }
144 
146  template<class Archive>
147  void serialize(Archive & ar, const unsigned int /*version*/) {
148  ar & mX;
149  ar & mY;
150  ar & mD;
151  ar & orthonormal;
152  ar & updateCount;
153  }
154 
155  inline bool isOrthonormal() const {
156  return orthonormal;
157  }
158 
159  inline bool isSquareRoot() const {
160  return m_squareRoot;
161  }
162 
164  inline void squareRoot(bool mode) {
165  if (mode == m_squareRoot)
166  return;
167 
168  if (mode)
169  mD.unaryExprInPlace(Eigen::internal::scalar_sqrt_op<Real>());
170  else
171  mD.unaryExprInPlace(Eigen::internal::scalar_abs2_op<Real>());
172 
173  m_squareRoot = mode;
174  }
175 
177  void orthonormalize() {
178  // Check if we have something to do
179  if (isOrthonormal()) return;
180 
181  // Orthonornalize (keeping in mind that our diagonal might be the square root)
182  if (m_squareRoot)
183  mD.unaryExprInPlace(Eigen::internal::scalar_abs2_op<Real>());
185  if (m_squareRoot)
186  mD.unaryExprInPlace(Eigen::internal::scalar_sqrt_op<Real>());
187 
188  // If we can linearly combine, use it to reduce the future amount of computation
189  if (fs->canLinearlyCombine()) {
190  mX = fs->linearCombination(mX, mY, 1);
191  mY = Eigen::Identity<Scalar>(mX->size(),mX->size());
192  }
193 
194  orthonormal = true;
195 
196  }
197 
199  void multiplyBy(Real alpha) {
200  if (m_squareRoot) {
201  if (alpha < 0)
202  KQP_THROW_EXCEPTION_F(out_of_bound_exception, "Cannot multiply a kernel operator by a negative value (%g)", %alpha);
203  alpha = std::sqrt(alpha);
204  }
205 
206  mD.unaryExprInPlace(Eigen::internal::scalar_multiple_op<Real>(alpha));
207  }
208 
209 
211  Real trace() const {
212  if (isOrthonormal())
213  return m_squareRoot ? mD.cwiseAbs2().sum() : mD.sum();
214 
215  if (m_squareRoot)
216  return fs->k(*mX, mY, mD).trace();
217 
218  return fs->k(*mX, mY, mD, *mX, mY, RealVector::Ones(mX->size())).trace();
219  }
220 
222  void traceNormalize() {
223  this->multiplyBy((Scalar)1 / this->trace());
224  }
225 
226 };
227 
228 
229 }
230 
231 #endif
232