Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
accumulator.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_ACCUMULATOR_BUILDER_H__
19 #define __KQP_ACCUMULATOR_BUILDER_H__
20 
21 #include <boost/static_assert.hpp>
22 
23 #include <kqp/alt_matrix.hpp>
24 #include <kqp/kernel_evd.hpp>
25 #include <kqp/evd_utils.hpp>
26 
27 namespace kqp{
28 
29 # include <kqp/define_header_logger.hpp>
30  DEFINE_KQP_HLOGGER("kqp.kevd.accumulator");
31 
41  template <typename Scalar, bool can_linearly_combine> class AccumulatorKernelEVD : public KernelEVD<Scalar> {};
42 
43  // Specialisation when we know how to combine linearly
44  template <typename Scalar> class AccumulatorKernelEVD<Scalar, true> : public KernelEVD<Scalar> {
45  public:
46  enum {
47  use_linear_combination = 1
48  };
49 
50  KQP_SCALAR_TYPEDEFS(Scalar);
51 
52  AccumulatorKernelEVD(const FSpaceCPtr &fs) : KernelEVD<Scalar>(fs), fMatrix(fs->newMatrix()) {
53  }
54 
55  virtual ~AccumulatorKernelEVD() {
56  }
57 
58 
59  virtual void _add(Real alpha, const FMatrixCPtr &mX, const ScalarAltMatrix &mA) override {
60  // Just add the vectors using linear combination
61  if(alpha < 0) {
62  KQP_HLOG_WARN_F("Downdating with accumulator (alpha=%g) [threshold=%g]", % -kqp::epsilon())
63  // Just ignore if small enough
64  if (-alpha < kqp::epsilon())
65  return;
66  KQP_THROW_EXCEPTION_F(not_implemented_exception, "Downdating with accumulator (alpha=%g)", % alpha);
67  }
68 
69  FMatrixPtr fm = this->getFSpace()->linearCombination(mX, mA, Eigen::internal::sqrt(alpha));
70  fMatrix->add(*fm);
71  }
72 
73  void reset() override {
75  fMatrix = this->getFSpace()->newMatrix();
76  }
77 
79  virtual Decomposition<Scalar> _getDecomposition() const override {
80 
81  const ScalarMatrix& gram = this->getFSpace()->k(fMatrix);
82  Eigen::SelfAdjointEigenSolver<ScalarMatrix> evd(gram.template selfadjointView<Eigen::Lower>());
83 
84  ScalarMatrix _mY;
85  RealVector _mD;
86  kqp::ThinEVD<ScalarMatrix>::run(evd, _mY, _mD);
87 
88  ScalarMatrix __mY;
89  __mY.noalias() = _mY * _mD.cwiseAbs().cwiseSqrt().cwiseInverse().asDiagonal();
90 
91  return Decomposition<Scalar>(this->getFSpace(), fMatrix, __mY, _mD, true, false);
92  }
93 
94  private:
96  FMatrixPtr fMatrix;
97  };
98 
99 
100 
101 
102  // Specialisation when we don't know how to combine linearly
103  template <typename Scalar> class AccumulatorKernelEVD<Scalar, false> : public KernelEVD<Scalar> {
104  public:
105 
106  KQP_SCALAR_TYPEDEFS(Scalar);
107 
108  AccumulatorKernelEVD(const FSpaceCPtr &fs) : KernelEVD<Scalar>(fs), fMatrix(fs->newMatrix()) {
109  offsets_X.push_back(0);
110  offsets_A.push_back(0);
111  }
112 
113  virtual ~AccumulatorKernelEVD() {}
114 
115  protected:
116  virtual void _add(Real alpha, const FMatrixCPtr &mX, const ScalarAltMatrix &mA) override {
117  if(alpha < 0)
118  KQP_THROW_EXCEPTION(not_implemented_exception, "Downdating with accumulator");
119 
120  // If there is nothing to add
121  if (mA.cols() == 0)
122  return;
123 
124  // Do a deep copy of mA
125  combination_matrices.push_back(mA);
126  alphas.push_back(Eigen::internal::sqrt(alpha));
127  fMatrix->add(mX);
128  offsets_X.push_back(offsets_X.back() + mX->size());
129  offsets_A.push_back(offsets_A.back() + mA.cols());
130  }
131 
132  void reset() {
134 
135  fMatrix = this->getFSpace()->newMatrix();
136  alphas.clear();
137  offsets_X.clear();
138  offsets_A.clear();
139  combination_matrices.clear();
140 
141  offsets_X.push_back(0);
142  offsets_A.push_back(0);
143  }
144 
145 
146  protected:
148  virtual Decomposition<Scalar> _getDecomposition() const override {
149  Decomposition<Scalar> d(this->getFSpace());
150  // Compute A^T X^T X A^T
151  // where A = diag(A_1 ... A_n) and X = (X_1 ... X_n)
152 
153  Index size = offsets_A.back();
154 
155  // Nothing to do
156  if (size == 0) {
157  d.mX = this->getFSpace()->newMatrix();
158  d.mY.resize(0,0);
159  d.mD.resize(0,1);
160  return d;
161  }
162 
163  ScalarMatrix gram_X = this->getFSpace()->k(fMatrix);
164  ScalarMatrix gram(size, size);
165 
166  for(size_t i = 0; i < combination_matrices.size(); i++) {
167  const ScalarAltMatrix &mAi = combination_matrices[i];
168  for(size_t j = 0; j <= i; j++) {
169  const ScalarAltMatrix &mAj = combination_matrices[j];
170  getBlock(gram, offsets_A, i, j)
171  = (mAi.adjoint() * ((Eigen::internal::conj(alphas[i]) * alphas[j]) * getBlock(gram_X, offsets_X, i, j))) * mAj;
172  }
173  }
174 
175 
176  // Direct EVD
177  ScalarMatrix _mY;
178  RealVector _mD;
179  Eigen::SelfAdjointEigenSolver<ScalarMatrix> evd(gram.template selfadjointView<Eigen::Lower>());
180  kqp::ThinEVD<ScalarMatrix>::run(evd, _mY, _mD);
181  d.mD.swap(_mD);
182 
183  // Y <- A * Y * D^-1/2 (use cwiseAbs to avoid problems with small negative values)
184  ScalarMatrix _mY2 = _mY * d.mD.cwiseAbs().cwiseSqrt().cwiseInverse().asDiagonal();
185  _mY = _mY2;
186  ScalarMatrix __mY(offsets_X.back(), _mY.cols());
187 
188  for(size_t i = 0; i < combination_matrices.size(); i++) {
189  const ScalarAltMatrix &mAi = combination_matrices[i];
190  __mY.block(offsets_X[i], 0, offsets_X[i+1]-offsets_X[i], __mY.cols()) = mAi * (alphas[i] * _mY.block(offsets_A[i], 0, offsets_A[i+1]-offsets_A[i], _mY.cols()));
191  }
192 
193  d.mY.swap(__mY);
194  d.mX = fMatrix;
195  return d;
196  }
197 
198  private:
199  static inline Eigen::Block<ScalarMatrix> getBlock(ScalarMatrix &m, const std::vector<Index> &offsets, size_t i, size_t j) {
200  return m.block(offsets[i], offsets[j], offsets[i+1] - offsets[i], offsets[j+1]-offsets[j]);
201  }
202 
204  FMatrixPtr fMatrix;
205 
207  std::vector<ScalarAltMatrix> combination_matrices;
208 
210  std::vector<Index> offsets_A;
211  std::vector<Index> offsets_X;
212 
213 
215  std::vector<Scalar> alphas;
216 
217  };
218 
219 }
220 
221 #endif