Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
probabilities.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_PROBABILITIES_H__
19 #define __KQP_PROBABILITIES_H__
20 
21 
22 #include <kqp/kqp.hpp>
23 #include <kqp/alt_matrix.hpp>
24 #include <kqp/kernel_evd.hpp>
25 #include <kqp/evd_utils.hpp>
26 
27 
28 namespace kqp {
29 
30  // Foward declarations
31  template <typename Scalar> class Density;
32  template <typename Scalar> class Event;
33 
34 
55  template <typename Scalar> class KernelOperator {
56  public:
57  KQP_SCALAR_TYPEDEFS(Scalar);
58 
64  KernelOperator(const KernelEVD<Scalar> & evd) : m_operator(evd.getDecomposition()) {
65  m_operator.squareRoot(true);
66  }
67 
75  m_operator.squareRoot(true);
76  }
77 
78 
82  KernelOperator(const FSpaceCPtr &fs, const FMatrixCPtr &mX, const ScalarAltMatrix &mY, const RealVector &mS, bool orthonormal)
83  : m_operator(fs, mX,mY,mS,orthonormal,false) {
84  m_operator.squareRoot(true);
85  }
86 
87 
95  void trim(Index newRank) {
96  this->orthonormalize();
97 
98  // Don't do anything
99  if (newRank <= getRank()) return;
100 
101  m_operator.mY.conservativeResize(m_operator.mY.rows(), newRank);
102  m_operator.mD.conservativeResize(newRank, 1);
103  }
104 
108  Index getRank() const {
109  if (isOrthonormal())
110  return m_operator.mD.size();
111  return -1;
112  }
113 
115  inline const FSpaceBase &fs() const {
116  return *m_operator.fs;
117  }
118 
120  inline const FMatrixBase &X() const {
121  return *m_operator.mX;
122  }
123 
124  inline const ScalarAltMatrix &Y() const {
125  return m_operator.mY;
126  }
127 
128  inline const RealAltVector &S() const {
129  return m_operator.mD;
130  }
131 
132  void orthonormalize() const {
133  const_cast<KernelOperator*>(this)->_orthonormalize();
134  }
135 
136  inline bool isOrthonormal() const {
137  return m_operator.orthonormal;
138  }
139 
143  Real squaredNorm() {
144  if (isOrthonormal())
145  return m_operator.mD.squaredNorm();
146 
147  return m_operator.fs->k(X(), Y(), S()).squaredNorm();
148  }
149 
156  ScalarMatrix inners(const KernelOperator<Scalar>& that) const {
157  return m_operator.fs->k(m_operator.mX, m_operator.mY, m_operator.mD,
158  that.m_operator.mX, that.m_operator.mY, that.m_operator.mD);
159  }
160 
165  FMatrixPtr matrix() const {
166  return m_operator.fs->linearCombination(X(), ScalarMatrix(Y() * S().asDiagonal()));
167  }
168 
169 
170  protected:
171 
174  m_operator.orthonormalize();
175  }
176 
179  friend class Density<Scalar>;
180  friend class Event<Scalar>;
181  };
182 
183 
184 
193  template <typename Scalar>
194  class Event : public KernelOperator<Scalar> {
195  bool useLinearCombination;
196  public:
197  KQP_SCALAR_TYPEDEFS(Scalar);
198 
204  Event(const KernelEVD<Scalar> &evd, bool fuzzy = false)
205  : KernelOperator<Scalar>(evd),
206  useLinearCombination(evd.getFSpace()->canLinearlyCombine()) {
207  if (!fuzzy) {
208  this->orthonormalize();
209  this->m_operator.mD = RealVector::Ones(this->Y().cols());
210  }
211  }
212 
214  Event(const FSpaceCPtr &fs, const FMatrixCPtr &mX, bool orthonormal)
215  : KernelOperator<Scalar>(fs, mX, Eigen::Identity<Scalar>(mX->size(),mX->size()), RealVector::Ones(mX->size()), orthonormal) {
216  }
217 
219  Event(const FSpaceCPtr &fs, const FMatrixCPtr &mX, const ScalarAltMatrix &mY, bool orthonormal)
220  : KernelOperator<Scalar>(fs, mX, mY, RealVector::Ones(mX->size()), orthonormal) {
221  }
222 
223  Event(const Decomposition<Scalar> & d, bool fuzzy) : KernelOperator<Scalar>(d) {
224  if (!fuzzy) {
225  this->orthonormalize();
226  this->m_operator.mD = RealVector::Ones(this->Y().cols());
227  }
228  }
229 
231  void setUseLinearCombination(bool b) {
232  this->useLinearCombination = b;
233  }
234 
235 
236 
237  friend class Density<Scalar>;
238 #ifndef SWIG
239  static Density<Scalar> projectWithLC(const Density<Scalar>& density, const Event<Scalar> &event) {
240  event.orthonormalize();
241  ScalarMatrix lc;
242  noalias(lc) = event.Y() * event.m_operator.k(density.m_operator);
243 
244  FMatrixPtr mX = event.m_operator.fs->linearCombination(event.X(), lc);
245  ScalarAltMatrix mY = Eigen::Identity<Scalar>(mX->size(),mX->size());
246  RealVector mS = RealVector::Ones(mY.cols());
247  return Density<Scalar>(event.m_operator.fs, mX, mY, mS, false);
248  }
249 
250  static Density<Scalar> projectOrthogonalWithLC(const Density<Scalar>& density, const Event<Scalar> &event) {
251  // We need the event to be in an orthonormal form
252  event.orthonormalize();
253 
254  Index n = event.S().rows();
255 
256  // FIXME
257  RealVector s = event.S();
258  ScalarAltMatrix e_mY(event.Y() * (RealVector::Ones(n) - (RealVector::Ones(n) - s.cwiseAbs2()).cwiseSqrt()).asDiagonal()
259  * event.m_operator.fs->k(event.X(), event.Y(), density.X(), density.Y()));
260 
261  FMatrixPtr mX = density.m_operator.fs->linearCombination(density.X(), density.Y(), 1., event.X(), e_mY, -1.);
262 
263  return Density<Scalar>(event.m_operator.fs, mX, Eigen::Identity<Scalar>(mX->size(),mX->size()), density.S(), false);
264  }
265 
266 
267  static Density<Scalar> project(const Density<Scalar>& density, const Event<Scalar> &event) {
268  event.orthonormalize();
269  FMatrixPtr mX = event.X().copy(); //event.m_operator.fs->newMatrix(event.X());
270  ScalarAltMatrix mY(event.Y() * event.m_operator.k(density.m_operator));
271  RealVector mS = RealVector::Ones(mY.cols());
272  return Density<Scalar>(event.m_operator.fs, mX, mY, mS, false);
273  }
274 
275  static Density<Scalar> projectOrthogonal(const Density<Scalar>& density, const Event<Scalar> &event) {
276  event.orthonormalize();
277 
278  // Concatenates the vectors
279  FMatrixPtr mX = density.X().copy(); //event.m_operator.fs->newMatrix(density.X());
280  mX->add(event.X());
281 
282  // Computes the new linear combination matrix
283  ScalarMatrix _mY(density.X().size() + event.X().size(), density.S().rows());
284  RealVector s = event.S();
285  Index n = s.rows();
286 
287  _mY.topRows(density.X().size()) = density.Y();
288  _mY.bottomRows(event.X().size()) = ((Scalar)-1)
289  * (event.Y() * (RealVector::Ones(n) - (RealVector::Ones(n) - s.cwiseAbs2()).cwiseSqrt()).asDiagonal() * event.m_operator.fs->k(event.X(), event.Y(), density.X(), density.Y()));
290 
291  ScalarAltMatrix mY;
292  mY.swap(_mY);
293 
294 
295  // Return
296  return Density<Scalar>(event.m_operator.fs, mX, mY, density.S(), false);
297  }
298 #endif
299 
304  Density<Scalar> project(const Density<Scalar>& density, bool orthogonal = false, bool normalize = true) {
305  // Orthogonal projection
306  if (orthogonal) {
307  Density<Scalar> r = useLinearCombination ? projectOrthogonalWithLC(density, *this) : projectOrthogonal(density, *this);
308  if (normalize)
309  r.normalize();
310  return r;
311  }
312 
313  // Normal projection
314  Density<Scalar> r = useLinearCombination ? projectWithLC(density, *this) : project(density, *this);
315  if (normalize)
316  r.normalize();
317  return r;
318  }
319  };
320 
326  template <typename Scalar> class Density: public KernelOperator<Scalar> {
327  public:
328  KQP_SCALAR_TYPEDEFS(Scalar);
329 
336  Density(const KernelEVD<Scalar> & evd) : KernelOperator<Scalar>(evd) {
337  }
338 
339  Density(const FSpaceCPtr &fs, const FMatrixCPtr &mX, const ScalarAltMatrix &mY, const RealVector &mD, bool orthonormal) : KernelOperator<Scalar>(fs, mX, mY, mD, orthonormal) {
340  }
341 
342  Density(const FSpaceCPtr &fs, const FMatrixCPtr &mX, bool orthonormal)
343  : KernelOperator<Scalar>(fs, mX, Eigen::Identity<Scalar>(mX->size(),mX->size()), RealVector::Ones(mX->size()), orthonormal) {
344  }
345 
346  Density(const Decomposition<Scalar> & d) : KernelOperator<Scalar>(d) {
347  }
348 
349 
351  void normalize() {
352  this->m_operator.traceNormalize();
353  }
354 
362  Real probability(const Event<Scalar>& event) const {
363  return this->m_operator.k(event.m_operator).squaredNorm();
364  }
365 
371  RealVector eigenProbabilities(const Event<Scalar>& event, bool noEigenValues) const {
372  this->orthonormalize();
373 
374  auto &other = event.m_operator;
375 
376  if (noEigenValues)
377  return this->m_operator.fs->k(
378  this->X(), this->Y(), RealVector::Ones(this->Y().cols()),
379  *other.mX, other.mY, other.mD)
380  .rowwise().squaredNorm();
381 
382  return this->m_operator.k(other).rowwise().squaredNorm();
383  }
384 
385 
387  Real entropy() const {
388  this->orthonormalize();
389  RealVector s = this->S();
390  return - (2 * s.array().log() * s.array().abs2()).sum();
391  }
392 
393  Real computeDivergence(const Density<Scalar> &tau) const {
394  return computeDivergence(tau, epsilon());
395  }
396 
405  Real computeDivergence(const Density<Scalar> &tau, Real epsilon) const {
406  const Density<Scalar> &rho = *this;
407 
408  // --- Requires orthonormal decompositions
409  rho.orthonormalize();
410  tau.orthonormalize();
411 
412  // --- Notation
413 
414  ScalarMatrix inners = this->m_operator.fs->k(rho.X(), tau.X());
415  noalias(inners) = rho.S().asDiagonal() * rho.Y().adjoint() * inners * tau.Y();
416 
417  // --- Compute tr(p log q)
418  Scalar plogq = 0;
419 
420 
421  // The background density span the subspace of rho
422  Index rank = rho.S().rows() + tau.S().rows();
423  Index dimension = rho.fs().dimension();
424  if (rank > dimension && dimension > 0) rank = dimension;
425 
426  Real alpha = 1. / (Real)(rank);
427  Real alpha_noise = epsilon * alpha;
428 
429  // Includes the smoothing probability if not too small
430  if (epsilon >= kqp::epsilon())
431  plogq = log(alpha_noise) * (1. - inners.squaredNorm());
432 
433 
434  // Main computation
435  RealVector mD(tau.S().rows());
436 
437  RealVector tau_S = tau.S();
438 
439  plogq -= (inners * (-((1 - epsilon) * tau_S.array().abs2() + alpha_noise).log()) .sqrt().matrix().asDiagonal()).squaredNorm();
440 
441 
442  // --- Compute tr(p log p)
443  return - plogq - this->entropy();
444 
445  }
446 
447 
448  };
449 
450 
451 
452 }
453 
454 #endif